CL reprojection sparse reconstruction method based on adaptive entropy weighting

Abstract

The invention discloses a CL reprojection sparse reconstruction method based on adaptive entropy weighting. The method comprises the following steps: 1, acquiring a group of original projection sequences P under an X-ray CL scanning system; 2, executing initialization operation, and reconstructing the projection sequence by adopting a SART three-dimensional iterative reconstruction algorithm under regularization constraint to obtain a reference image f * as a reference quantity; 3, performing re-projection operation on f * to obtain a projection sequence P'of reconstructed voxels; 4, subtracting P from P 'to obtain a differential projection sequence delta P with a relatively high sparse degree; 5, reconstructing the sparse delta P by using the l1 norm minimization to obtain a sparse reconstructed image; and 6, fusing the f * and the sparse reconstructed image by using the entropy coefficient of the corresponding layer of the f * and the sparse reconstructed image. According to the method, the interlayer aliasing problem caused by a CL scanning mode under complete angle reconstruction or incomplete angle reconstruction can be greatly reduced under the condition that the number of iterations is small, and therefore the reconstruction precision is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of CL image reconstruction, and in particular to a CL reprojection sparse reconstruction method based on adaptive entropy value weighting. Background Art

[0002] In industrial nondestructive testing, the CT (Computed Tomography) method cannot be effectively applied to the 3D reconstruction of flat objects, so the CL (Computed Laminography) method is usually used to reconstruct flat objects in 3D. However, the projection data collected by the CL scanning method is not as complete as that of the CT scanning method. The use of traditional reconstruction methods such as (SART, SART_TV, etc.) will cause certain artifacts in the 3D reconstruction results, which are specifically manifested as (X s O s Y s plane) has a good effect, but there is aliasing between slices (X s O s Z s and Y s O s Z s plane). The reconstruction reference difference algorithm (RRD algorithm) was originally proposed for the reconstruction of incomplete CT angles. The algorithm uses the reconstruction result of differential sparse projection and the reference image reconstructed by the ART algorithm to sum, thereby improving the overall accuracy of the reconstructed image. However, for CL reconstruction with good effects within slices and serious aliasing between slices, simple summation cannot effectively reduce the inter-layer aliasing phenomenon and can only improve the reconstruction accuracy within the slice to a certain extent. Other reconstruction algorithms that add various constraints in the iterative reconstruction process, such as (TV regularization constraints, RTV relative total variation), usually require the definition of multiple hyperparameters to adjust the reconstruction quality. The algorithm has poor robustness and requires a high number of iterations, resulting in low reconstruction efficiency. Summary of the invention

[0003] In order to solve the deficiencies of the above-mentioned prior art, the present invention proposes a CL reprojection sparse reconstruction method based on adaptive entropy value weighting, so as to effectively eliminate the inter-layer aliasing phenomenon generated in the CL reconstruction process while also improving the image reconstruction efficiency.

[0004] In order to achieve the above-mentioned purpose, the present invention adopts the following technical solutions:

[0005] The CL reprojection sparse reconstruction method based on adaptive entropy value weighting of the present invention is characterized in that it is applied to the three-dimensional reconstruction process of the CL system and is performed according to the following steps:

[0006] Step 1: Take the rotation center of the plate as the origin , with the horizontal direction of the plate plane as X s Axis direction, the longitudinal direction of the plate plane is Y s Axis direction, perpendicular to plane X s O s Y s The direction is Z s Axis direction, thereby establishing the object coordinate system O s X s Y s Z s , where Z s is the object's rotation axis;

[0007] Define the object coordinate system O s X s Y s Z s Plane X in s O s Y s is the direction within the slice, plane X s O s Z s With plane Y s O s Z s is the inter-slice direction.

[0008] Step 2: Use the CL system to collect the projection data set P of the image to be reconstructed at θ angles ,in, Represents the projection data of the image to be reconstructed at the kth angle, let The corresponding system projection matrix is ,and , Represents a system projection matrix collection.

[0009] Step 3: Define the current number of external iterations as n, and initialize n=0, and define the total number of external iterations as ; Define the iteration factor as ,and ;

[0010] Define the current number of inner iterations of the TV constraint process as , the total number of internal iterations of the TV constraint process is ; Define the TV constraint factor as ,and ;

[0011] Define and initialize the image to be reconstructed under the nth outer iteration .

[0012] Step 4: Use equation (1) to obtain the image to be reconstructed under the n+1th outer iteration :

[0013] (1)

[0014] In formula (1), and Respectively The sum of the elements in the i-th row and the sum of the elements in the j-th column; T represents transpose.

[0015] Step 5: The pixels with grayscale values ​​less than 0 are set to 0, thus obtaining the image to be reconstructed after the n+1th outer iteration. .

[0016] Step 6: Initialization ; Define the n+1th outer iteration The image to be reconstructed in the inner iteration is recorded as , and initialize .

[0017] Step 7: Use formula (2) to Make corrections and get the n+1th outer iteration The image to be reconstructed under the inner iteration :

[0018] (2)

[0019] In formula (2), for and The difference, express The gradient image.

[0020] Step 8: Assign to Then return to step 7 and execute sequentially until So far, we get the n+1th outer iteration The image to be reconstructed under the inner iteration , and recorded as the updated image to be reconstructed under the n+1th outer iteration And assign it to .

[0021] Step 9: After assigning n to n+1, return to step 4 and execute sequentially until n>N, thereby obtaining the image to be reconstructed under the Nth outer iteration. , and recorded as the reference image .

[0022] Step 10: According to the geometric parameters of the CL system, use equation (3) to obtain the angle θ The projection data set ;

[0023] (3)

[0024] In formula (3), Indicates the kth angle The projection data.

[0025] Step 11: Use formula (4) to obtain the differential projection set ;

[0026] (4)

[0027] In formula (4), Represents the differential projection data at the kth angle.

[0028] Step 12: Use formula (5) to build a sparse reconstruction model And solve it to get the sparse reconstructed image ;

[0029] (5)

[0030] In formula (5), express norm, is the optimization termination threshold.

[0031] Step 13: Define the reference image Reconstruct the image with coefficients Respectively in Z s The image in the M-layer slice in the axial direction is recorded as and ,in, express In Z s The image in the mth slice in the axial direction, express In Z s The image of the mth slice in the axial direction.

[0032] Step 14: Calculate using formula (6) In Z s The mth layer image in the axis direction The entropy value as well as In Z s The mth layer image in the axis direction The entropy value :

[0033] (6)

[0034] In formula (7), and Slice images and The maximum gray value in Indicates that the pixel with gray value a is Probability of occurrence; Indicates that the pixel with gray value a is Probability of occurrence.

[0035] Step 15: Use equation (7) to get Z s The mth layer reconstructed image in the axial direction , thus obtaining the object to be measured at Z s M-layer reconstructed image set in the axial direction .

[0036] (7)

[0037] The electronic device of the present invention includes a memory and a processor, and is characterized in that the memory is used to store a program that supports the processor to execute the adaptive entropy value weighted CL reprojection sparse reconstruction method, and the processor is configured to execute the program stored in the memory.

[0038] The present invention provides a computer-readable storage medium, wherein a computer program is stored on the computer-readable storage medium, and the computer program, when executed by a processor, executes the steps of the adaptive entropy value weighted CL reprojection sparse reconstruction method.

[0039] Compared with the prior art, the present invention is beneficial in that:

[0040] 1. The present invention adopts CL to reconstruct image features and uses the entropy values ​​of slice images of different layers as adaptive weighting coefficients to fuse the reference image with the differential vector, so that it can eliminate the inter-layer aliasing phenomenon caused by the CL scanning method to a large extent.

[0041] 2. The present invention uses the reconstruction result of the SART_TV algorithm as a reference image to increase the sparsity of the differential projection, thereby further improving the final reconstruction quality.

[0042] 3. Iterative 3D reconstruction methods that include various constraints usually include the definition of many hyperparameters. The method adopted in the present invention does not require additional hyperparameters to be defined, so it has better robustness.

[0043] 4. The reconstruction method proposed in the present invention no longer requires multiple iterations to achieve a better reconstruction effect, which effectively improves the reconstruction speed. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] Figure 1A geometric scanning model diagram of a CL system in an embodiment of the present invention;

[0045] Figure 2 This is a flowchart of the algorithm implementation in the embodiment of the present invention;

[0046] Figure 3 A model diagram of theoretical numerical simulation in an embodiment of the present invention;

[0047] Figure 4 A slice diagram of a reference image and a sparsely reconstructed image in an embodiment of the present invention;

[0048] Figure 5 This is a comparison chart of the reconstruction effects of the SART algorithm (SART_TV) under regularization constraints under complete acquisition angle reconstruction in an embodiment of the present invention, the RRD algorithm, and the algorithm proposed by the present invention;

[0049] Figure 6 This is a comparison chart of the reconstruction effects of the SART algorithm (SART_TV) under regularization constraints under incomplete acquisition angle reconstruction in an embodiment of the present invention, the RRD algorithm, and the algorithm proposed in the present invention. DETAILED DESCRIPTION

[0050] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0051] In this embodiment, in order to reduce the inter-layer aliasing problem caused by the CL scanning method in the 3D reconstruction process, a CL reprojection sparse reconstruction method based on adaptive entropy value weighting is proposed, which uses the SART_TV algorithm to replace the ART algorithm in the RRD algorithm to reconstruct a reference image that is closer to the original reconstructed object, improve the sparsity of the differential projection, and use The optimization criterion of the minimum norm is used to reconstruct the differential projection to obtain a sparse reconstructed image. Finally, the entropy value of the slice image of different layers in the Zs axis direction of the reconstructed image is used as a weighting coefficient to fuse the reference image with the sparse reconstructed image, thereby reducing the degree of inter-layer aliasing in the CL reconstruction process. Specifically, the method includes the following steps:

[0052] Step 1: Take the rotation center of the plate as the origin , with the horizontal direction of the plate plane as X s Axis direction, the longitudinal direction of the plate plane is Y s Axis direction, perpendicular to plane X s O s Y s The direction is Z s Axis direction, thereby establishing the object coordinate system O s X s Y s Z s , where Zs is the object's rotation axis;

[0053] Define the object coordinate system O s X s Y s Z s Plane X in s O s Y s is the direction within the slice, plane X s O s Z s With plane Y s O s Z s is the inter-slice direction;

[0054] The center of the detector is taken as the origin , with the horizontal direction of the detector plane as the X-axis direction, the vertical direction of the detector plane as the Z-axis direction, and the direction perpendicular to the detector plane as the Y direction, to establish the world coordinate system OXYZ;

[0055] Define the object coordinate system O s X s Y s Z s The object rotation axis Z s The angle between the axis α and the Z axis of the world coordinate system OXYZ is α, and α is the inclination angle of the CL system. In this embodiment, the geometry of the CL system is as follows: Figure 1 As shown, the geometric parameters of the system are as follows: the distance from the ray source to the object is 1131.37 mm, the distance from the ray source to the center of the detector is 2404.16 mm, the sample size is 323 × 378 × 102 mm, and the voxel size is 1 mm. Figure 3 As shown, the size of the detector is 500mm×500mm, the pixel size is 1mm; the tilt angle of the CL system is 45°; under the geometric relationship of this system, according to Figure 2 The steps shown implement the method proposed by the present invention.

[0056] Step 2: Use the CL system to collect the projection data set P of the object to be measured at θ angles ,in, Represents the projection data of the object to be measured at the kth angle, let The corresponding system projection matrix is ,and , represents the projection matrix set, Figure 3 The PCB original model in the embodiment and its Figure 1 A collection of projection images under the CL system shown.

[0057] Step 3: Define the current number of external iterations as n, and initialize n=0, and define the total number of external iterations as ; Define the iteration factor as ,and In this embodiment, ;

[0058] Define the current number of inner iterations of the TV constraint process as , the total number of internal iterations of the TV constraint process is ; Define the TV constraint factor as ,and In this embodiment, ;

[0059] Define and initialize the image to be reconstructed under the nth outer iteration .

[0060] Step 4: Use equation (1) to obtain the image to be reconstructed under the n+1th outer iteration :

[0061] (1)

[0062] In formula (1), and Respectively The sum of the elements in the i-th row and the sum of the elements in the j-th column; T represents transpose;

[0063] The SART algorithm updates the pixel value of the object to be measured by calculating all rays passing through the pixel at the same projection angle. It is worth noting that in formula (1) and It appears in the denominator, so its value must not be 0 in actual implementation.

[0064] Step 5: The pixels with grayscale values ​​less than 0 are set to 0, thus obtaining the image to be reconstructed after the n+1th outer iteration. , that is, using formula (2) Perform non-negativity constraints.

[0065] (2)

[0066] Step 6: Initialization ; Define the n+1th outer iteration The image to be reconstructed in the inner iteration is recorded as , and initialize .

[0067] Step 7: Use formula (3) to Make corrections and get the n+1th outer iteration The image to be reconstructed under the inner iteration :

[0068] (3)

[0069] In formula (3), for and The difference, express The calculation method of the gradient image is shown in formula (4) and formula (5);

[0070] (4)

[0071] (5)

[0072] In formula (4) and formula (5), , and Respectively Index values ​​in the three axis directions, Usually a very small value is taken. In the embodiment, .

[0073] Step 8: Assign to Then return to step 7 and execute sequentially until So far, we get the n+1th outer iteration The image to be reconstructed under the inner iteration The updated image to be reconstructed under the n+1th outer iteration is recorded as And assign it to .

[0074] Step 9: After assigning n to n+1, return to step 4 and execute sequentially until n>N, thereby obtaining the image to be reconstructed under the Nth outer iteration. , and recorded as the reference image .

[0075] Step 10: According to the geometric parameters of the CL system, use equation (6) to obtain the angle θ The projection data set , since the reference image is very close to the original object, so and It’s also very close;

[0076] (6)

[0077] In formula (6), Indicates the kth angle Projection data, initial projection process and reference image The projection process must keep the projection matrix W consistent.

[0078] Step 11: Due to and Very close, according to the compressed sensing theory, in this embodiment, a differential projection vector with a very high degree of sparseness can be constructed, and then the differential projection vector can be reconstructed using the reconstruction algorithm, thereby obtaining a reconstruction result with less inter-layer aliasing. Use formula (7) to obtain the differential projection set ;

[0079] (7)

[0080] In formula (7), Represents the differential projection data at the kth angle.

[0081] Step 12: Use formula (8) to construct a sparse reconstruction model And solve it using The optimal sparse reconstructed image is obtained by the minimum norm method ;

[0082] (8)

[0083] In formula (8), express norm, is the optimization termination threshold. In the embodiment, , The norm refers to the sum of the absolute values ​​of all elements in a vector. Its advantage is that it can produce a sparse solution, that is, most of the elements in the vector are set to 0. The present invention uses the fastest gradient descent method to solve the optimization problem of formula (8). Figure 4 For this embodiment, the reference image and sparsely reconstructed images In Z s Comparison of the 46th and 50th slice images in the axial direction shows that the reference image The slice image contains more detailed information, but there is inter-layer aliasing, while the sparsely reconstructed image The slice image mainly retains the edge information and has no inter-layer aliasing phenomenon.

[0084] Step 13: Define the reference image and the coefficient reconstruction vectors are respectively in Z s The image in the M-layer slice in the axial direction is recorded as and ,in, express In Zs The image in the mth slice in the axial direction, express In Z s The image of the mth slice in the axial direction.

[0085] Step 14: Calculate using formula (9) In Z s The mth layer image in the axis direction The entropy value as well as In Z s The mth layer image in the axis direction The entropy value :

[0086] (9)

[0087] In formula (9), and Slice images and The maximum gray value in Indicates that the pixel with gray value a is Probability of occurrence; Indicates that the pixel with gray value a is The probability of occurrence, and Represents the reference images and sparsely reconstructed images In Z s The image entropy value of the mth layer in the axial direction.

[0088] Step 15: From Figure 4 As can be seen from the results, when the reconstructed image is at Z s When the slices in the axial direction contain more detail information, the image entropy value of the reference image is larger. At this time, there are fewer slice images corresponding to the reference image in the sparsely reconstructed image, which are mainly edge information. Therefore, the slice image entropy value of the reference image is used as a weighting coefficient to multiply the sparsely reconstructed image, and the sum of the two can enhance the edge features within the slice, thereby reducing aliasing. However, in actual situations, if the entropy value of the slice image of the reference image is small, it may cause the slice image pixel values ​​of the coefficient reconstructed image to be all 0, which in turn causes the entropy value of the slice image to be 0, causing the fusion result to lose some detail information. Therefore The optimal reconstruction with the minimum norm also requires adjusting the number of iterative reconstruction algorithms to avoid losing too much detail information. This needs to be determined according to the different structures of the object. s The entropy value of each layer of the image in the axial direction is the reference image and sparsely reconstructed images The fusion method is to use the reference image The entropy value of the m-layer image is used as the sparse reconstructed image The weights of the m-layer image are reconstructed using sparse The entropy value of the m-layer image is used as the reference image The weights of the m layers of images are summed. That is, using formula (10) we get Z s The mth layer reconstructed image in the axial direction , thus obtaining the object to be measured at Z s M-layer reconstructed image set in the axial direction ;

[0089] (10)

[0090] In this embodiment, for Figure 3 The PCB model shown in the figure is projected at full angles and incomplete angles respectively. The projection data acquisition angle at full angles is 0-359°( ), a total of 360 sets of projection data, the collection angle of the projection data under incomplete angles is 0-149°( ), a total of 150 angles. Steps 1 to 15 of the method proposed in the present invention are performed on the two sets of projection data respectively to obtain a reconstructed image and compare it with the SART_TV algorithm and the RRD algorithm. Figure 5 This is the comparison result of the reconstructed image at the full angle. Figure 6 The comparison results of the reconstructed images under incomplete angles are shown in Figure 2. Compared with the results of the original PCB model image, it can be seen that the method proposed by the present invention can largely eliminate the inter-layer aliasing phenomenon generated by the reconstruction process of the CL system under complete and incomplete angles. This effect is more obvious in the background layer slice image such as the 50th layer.

[0091] In this embodiment, an electronic device includes a memory and a processor, wherein the memory is used to store a program that supports the processor to execute the above method, and the processor is configured to execute the program stored in the memory.

[0092] In this embodiment, a computer-readable storage medium stores a computer program on the computer-readable storage medium, and the computer program executes the steps of the above method when executed by a processor.

Claims

1. A CL reprojection sparse reconstruction method based on adaptive entropy value weighting, characterized in that: It is applied in the 3D reconstruction process of the CL system and is carried out as follows: Step 1: Take the rotation center of the plate as the origin , with the horizontal direction of the plate plane as X s Axis direction, the longitudinal direction of the plate plane is Y s Axis direction, perpendicular to plane X s O s Y s The direction is Z s Axis direction, thereby establishing the object coordinate system O s X s Y s Z s , where Z s is the object's rotation axis; Define the object coordinate system O s X s Y s Z s Plane X in s O s Y s is the direction within the slice, plane X s O s Z s With plane Y s O s Z s is the inter-slice direction; Step 2: Use the CL system to collect the projection data set P of the image to be reconstructed at θ angles ,in, Represents the projection data of the image to be reconstructed at the kth angle, let The corresponding system projection matrix is ,and , Represents the system projection matrix set; Step 3: Define the current number of external iterations as n, and initialize n=0, and define the total number of external iterations as ; Define the iteration factor as ,and ; Define the current number of inner iterations of the TV constraint process as , the total number of internal iterations of the TV constraint process is ; Define the TV constraint factor as ,and ; Define and initialize the image to be reconstructed under the nth outer iteration ; Step 4: Use equation (1) to obtain the image to be reconstructed under the n+1th outer iteration : (1) In formula (1), and Respectively The sum of the elements in the i-th row and the sum of the elements in the j-th column; T represents transpose; Step 5: The pixels with grayscale values ​​less than 0 are set to 0, thus obtaining the image to be reconstructed after the n+1th outer iteration. ; Step 6: Initialization ; Define the n+1th outer iteration The image to be reconstructed in the inner iteration is recorded as , and initialize ; Step 7: Use formula (2) to Make corrections and get the n+1th outer iteration The image to be reconstructed under the inner iteration : (2) In formula (2), for and The difference, express The gradient image of Step 8: Assign to Then return to step 7 and execute sequentially until So far, we get the n+1th outer iteration The image to be reconstructed under the inner iteration , and recorded as the updated image to be reconstructed under the n+1th outer iteration And assign it to ; Step 9: After assigning n to n+1, return to step 4 and execute sequentially until n>N, thereby obtaining the image to be reconstructed under the Nth outer iteration. , and recorded as the reference image ; Step 10: According to the geometric parameters of the CL system, use equation (3) to obtain the angle θ The projection data set ; (3) In formula (3), Indicates the kth angle Projection data of Step 11: Use formula (4) to obtain the differential projection set ; (4) In formula (4), represents the differential projection data at the kth angle; Step 12: Use formula (5) to build a sparse reconstruction model And solve it to get the sparse reconstructed image ; (5) In formula (5), express norm, is the optimization termination threshold; Step 13: Define the reference image Reconstruct the image with coefficients Respectively in Z s The image in the M-layer slice in the axial direction is recorded as and ,in, express In Z s The image in the mth slice in the axial direction, express In Z s The image in the mth slice in the axial direction; Step 14: Calculate using formula (6) In Z s The mth layer image in the axis direction The entropy value as well as In Z s The mth layer image in the axis direction The entropy value : (6) In formula (7), and Slice images and The maximum gray value in Indicates that the pixel with gray value a is Probability of occurrence; Indicates that the pixel with gray value a is Probability of occurrence; Step 15: Use equation (7) to get Z s The mth layer reconstructed image in the axial direction , thus obtaining the object to be measured at Z s M-layer reconstructed image set in the axial direction ; (7)。 2. An electronic device, comprising a memory and a processor, characterized in that: The memory is used to store a program that supports the processor to execute the adaptive entropy value weighted CL reprojection sparse reconstruction method as described in claim 1, and the processor is configured to execute the program stored in the memory.

3. A computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the steps of the adaptive entropy value weighted CL reprojection sparse reconstruction method according to claim 1 are executed.

Images