An Adaptive Recognition Pulse-Type Terahertz CT Data Preprocessing Method

Abstract

The present invention relates to a method for preprocessing pulsed terahertz CT data with adaptive recognition. At each fixed coordinate position, multiple time-domain pulse signals are collected in a round-trip scanning mode, and the time-domain pulse signals obtained from multiple repeated measurements at the same coordinate position are averaged; peak-to-peak extraction based on an adaptive threshold is performed; the angular differences between each coordinate point and its adjacent points are detected, and the coordinate points meeting the interpolation conditions are selected. A linear interpolation method is applied between the coordinate points meeting the interpolation conditions to generate uniformly distributed interpolation coordinates, thereby realizing resampling of the signals; the peak-to-peak values of the signals are normalized, and the signals are converted from peak-to-peak values to absorption coefficients; the even rows are flipped; automatic calibration is achieved through a data matching method to adjust the rotation center of the sinogram. The adaptive threshold method is used to accurately extract the peak-to-peak values of the signals, ensuring a reliable characterization of the energy attenuation characteristics. A linear interpolation algorithm is introduced to unify the time and space scales and eliminate acquisition errors.

Description

Technical Field

[0001] The present invention belongs to the technical field of terahertz CT data processing, and particularly relates to a pulsed terahertz CT data preprocessing method for adaptive recognition. Background Art

[0002] X-Ray Computed Tomography (X-Ray CT) was initially developed as a medical imaging technology for generating tomographic images or "slices" of a specific area of a patient. These cross-sectional images are of great value in multiple medical fields, especially in diagnosis and treatment. X-Ray CT can reconstruct a true three-dimensional image of the interior of an object through a series of two-dimensional X-ray images taken from a single rotation axis (the emitter and sensor rotate around the patient). Nowadays, X-Ray CT is also widely used in the industrial field and has become an important tool for quality control and defect location. However, due to the certain biological hazards of X-ray radiation, its scope of application is significantly limited.

[0003] In 2002, Fergusson et al. introduced Terahertz Computed Tomography (Terahertz CT). Compared with X-ray imaging, Terahertz CT shows unique advantages in measuring amplitude and spectral phase information. Terahertz CT can identify or locate different substances in a non-destructive manner and has low radiation hazards, which makes it have great application potential in material analysis and biomedical imaging. Recently, Brahm et al. used Terahertz CT technology to achieve spectral analysis of materials and successfully distinguished lactose and glucose embedded in a polystyrene block through significant absorption characteristic lines in the terahertz spectrum. However, Terahertz CT technology still faces some challenges, especially the limitation of the absorption phenomenon, which restricts the thickness of the imaged sample and has become the main bottleneck for the transformation of Terahertz CT technology from the laboratory to practical applications.

[0004] In addition, Terahertz CT is also affected by equipment and environmental noise and unstable factors during the data acquisition process, which may reduce the accuracy of the imaging results. Therefore, in order to overcome the limitations brought by these challenges and noise interference, the preprocessing stage of terahertz (THz) CT data is particularly important in the terahertz CT system. Data preprocessing provides more reliable and accurate input data for image reconstruction by optimizing signals, correcting errors, and suppressing noise.

[0005] Existing THz data preprocessing methods are usually designed temporarily for specific tasks and lack a unified, systematic, and standardized framework. This design pattern leads to different researchers adopting diverse preprocessing steps according to their respective experimental requirements, resulting in diversity and inconsistency in methods, thereby weakening the comparability of data. In addition, most existing THz data preprocessing technologies rely heavily on manual experience. For example, in key steps such as noise removal, signal calibration, and offset correction, a large amount of manual intervention is often required. This not only reduces the processing efficiency but also easily introduces human errors, resulting in insufficient accuracy and reliability of the processing results.

[0006] Therefore, how to improve the existing THz data preprocessing methods, reduce manual intervention, and enhance the efficiency and accuracy of core steps such as noise removal, signal calibration, and offset correction, while improving data quality and ensuring the accuracy and stability of the results during the image reconstruction process, is a technical problem that urgently needs to be solved at present. Summary of the Invention

[0007] The purpose of the present invention is to provide an adaptive recognition pulsed terahertz CT data preprocessing method, which improves the existing THz data preprocessing method, reduces manual intervention, and enhances the efficiency and accuracy of core steps such as noise removal, signal calibration, and offset correction, while improving data quality and ensuring the accuracy and stability of the results during the image reconstruction process.

[0008] To solve the above technical problems, the technical solutions adopted by the present invention are as follows:

[0009] An adaptive recognition pulsed terahertz CT data preprocessing method includes the following steps:

[0010] S1: At each fixed coordinate position, multiple time-domain pulse signals with slight differences are collected through a round-trip scanning mode, and the time-domain pulse signals of multiple repeated measurements at the same coordinate position are averaged.

[0011] S2: Extract the peak-to-peak value of the signal based on the peak-to-peak value extraction method with an adaptive threshold.

[0012] S3: Detect the angular difference between each coordinate point and its adjacent points, select the coordinate points that meet the interpolation conditions, and apply the linear interpolation method between the coordinate points that meet the interpolation conditions to generate uniformly distributed interpolation coordinates to achieve resampling of the signal.

[0013] S4: Normalize the peak-to-peak value of the collected signal, convert the signal from the peak-to-peak value to the absorption coefficient, and calculate and correct the absorption coefficient of the signal.

[0014] S5: Flip the data corresponding to the even rows, that is, the data of the reverse scan.

[0015] S6: Automatic calibration is achieved through data matching method to adjust the rotation center of the sine graph.

[0016] Preferably, the specific formula for performing averaging processing on the time domain pulse signals repeatedly measured at the same coordinate position in step S1 is as follows:

[0017]

[0018] in, is the averaged signal, k is the total number of averaged signals, i is the number of a single signal, S i For a single signal.

[0019] Preferably, the specific process of performing peak-to-peak value extraction of the signal based on the peak-to-peak value extraction method based on the adaptive threshold in step S2 is as follows:

[0020] S21: Extract each time domain signal S i (t), and calculate each time domain signal S i (t), and calculate the peak-to-peak value of the signal. i The specific formulas for the maximum and minimum values of (t) are as follows:

[0021] max i =max(S i (t));

[0022] min i =min(S i (t));

[0023] Among them, max(·) is the maximum value function, and min(·) is the minimum value function;

[0024] S22: Introduce a vertical amplitude threshold to limit the peak amplitude, and combine it with a horizontal distance threshold to constrain the time interval between peaks;

[0025] S23: For signals exceeding the threshold limit, the peak-to-peak value extraction process is optimized by using the local extreme value search method, which optimizes the peak-to-peak value calculation by searching for local maximum and local minimum values within a window range;

[0026] S24: Define the search window max and Window min The window size is determined by the horizontal distance threshold x Decision, in the window Window max Find the local minimum value within min , and in the window Window minSearch for local maximum values max , and obtain the final peak-to-peak value.

[0027] Preferably, the specific process of introducing a vertical amplitude threshold to limit the peak amplitude and combining a horizontal distance threshold to constrain the time interval between peaks in step S22 is as follows:

[0028] When the absolute value exceeds the vertical amplitude threshold, i.e., |max i | > threshold y , or |min i | > threshold y , then continue to calculate the peak-to-peak value. If it does not exceed, it is regarded as the background noise of the signal; the horizontal distance threshold is constrained to threshold x ;

[0029] If the index difference between the maximum value t max and the minimum value t min is less than the horizontal distance threshold threshold x , i.e., |t max -t min | < threshold x , then it is determined that the maximum and minimum values come from a pulse signal under the horizontal distance threshold threshold x , and directly calculate the peak-to-peak value, Signal_FF = max i -min i .

[0030] Preferably, the specific formulas for defining the search windows Window max and Window min in step S24 are as follows:

[0031] Window max = [tma x -(threshold x -1), t max +(threshold x -1)];

[0032] Window min = [t min -(threshold x -1), t min +(threshold x -1)];

[0033] Search for local minimum values min , and search for local maximum values within the window Window min ​max The specific formula for calculating the final peak-to-peak value is as follows:

[0034] Signal_FF = max(|max i | + |local min |, |min i | + |local max |);

[0035] Where, Window max is the maximum value of the search window on the horizontal time axis, Window min is the minimum value of the search window on the horizontal time axis, Signal_FF is the peak-to-peak value of the time-domain pulse signal, max i is the maximum value of the vertical coordinate within the search window, min i is the maximum value of the vertical coordinate within the search window, t max is the horizontal coordinate corresponding to max i and t min is the horizontal coordinate corresponding to min i .

[0036] Preferably, in step S3, a linear interpolation method is applied between the coordinate points that meet the interpolation conditions to generate uniformly distributed interpolation coordinates. The specific process is as follows:

[0037] S31: For the uneven coordinates x and the corresponding signals y at each angle, the maximum travel in the y direction is y max . First, generate equally spaced horizontal coordinates x p , and the horizontal coordinates are equally spaced values in the interval [-y max , y max . The number of points is controlled by the resolution r, that is

[0038]

[0039] Where, is the interval between the horizontal coordinates, and is the value of the horizontal coordinate to be interpolated, k is the number of coordinates corresponding to the coordinates to be interpolated, and Δx is the equal interval distance of the coordinates to be interpolated;

[0040] S32: Use linear interpolation to calculate the corresponding signal value y p , and the specific formula is as follows:

[0041]

[0042] Where, is the difference function, and the corresponding signal value is calculated through the input horizontal coordinate ​ interpolation is the interpolation function, and linear is the linear interpolation function.

[0043] Preferably, the specific process of calculating and correcting the absorption coefficient of the signal in step S4 is as follows:

[0044] S41: Normalize the peak-to-peak value of the collected signal to obtain the transmittance. The specific calculation formula is as follows:

[0045]

[0046] S42: Calculate the absorbance of the signal. The specific formula is as follows:

[0047]

[0048] Among them, A i is the absorbance of the i-th signal point, Si is the peak-to-peak value of the i-th signal, s i ′ is the normalized peak-to-peak value;

[0049] S43: Define the correction factor ∈, when A i When it is less than ∈, force its value to be no less than ∈, that is, A i =max(A i ,∈).

[0050] Preferably, the specific process of performing the flipping process on the even-numbered rows, i.e., the data corresponding to the reverse scan, in step S5 is as follows:

[0051] By traversing the sinogram matrix, an element-by-element inversion operation is performed on each row of data starting from the second row of data, that is, the data with index 1, so that the direction of the reverse scan data is consistent with the forward scan data of the previous row.

[0052] Preferably, in step S6, automatic calibration is achieved by a data matching method, and the specific process of adjusting the rotation center of the sinusoidal graph is as follows:

[0053] S61: Let the data of the given sinusoidal graph be S(x, θ), where x represents the position of the detector and θ represents the projection angle; for each row of data, use the differential filter D = [1, -1] to calculate the difference of each row. The specific formula is as follows:

[0054] ΔS(x)=|D*S(x,θ)|;

[0055] Among them, * is the convolution operation, ΔS(x) represents the rate of change of two adjacent data;

[0056] S62: Normalize the intensity value of each row of data to extract the boundary points. The normalization of the i-th row of data can be expressed as:

[0057]

[0058] wherein, ∈ is a specified constant, and S(x, θ i ) is the projection data at the i-th angle, and S min (x, θ i ) is the minimum value of the projection data at the i-th angle, and S max (x, θ i ) is the maximum value of the projection data at the i-th angle;

[0059] S63: By calculating the rate of change ΔS(x) and the normalized intensity S norm (x, θ i ), for each row of projection data, extract the points where the data first and last change as the left boundary point x0 and the right boundary point x1;

[0060] S64: After extracting the left and right boundary points from each row, in order to minimize the deviation of the projection rotation center while retaining the position information of the measurement object, sort the left boundary point x0 and the right boundary point x1 in ascending and descending order respectively, and take the mean of the first boundary points as the optimal boundary points for matching and

[0061] S65: Calculate the offset Δ between the left and right boundary points. By calculating the offset Δ, fill the background data on both sides of the sinogram for compensation, and calibrate the intensity of the specified boundary points. The formula for calculating the offset Δ is as follows:

[0062]

[0063] Preferably, the specific formulas for extracting the points where the data first and last change as the left boundary point x0 and the right boundary point x1 in step S63 are as follows:

[0064] x0 = argmin x (ΔS(x) > Q3(ΔS(x)) ∩ S norm (x, θ i ));

[0065] x1 = argmax x (ΔS(x) > Q3(ΔS(x)) ∩ S norm (x, θ i ));

[0066] wherein, argmin x returns the position where x is the smallest and satisfies the condition, and argmax xTo return the position where the condition is satisfied and x is the largest, ΔS(x) is the change rate between two adjacent data, Q3 is the third quartile of ΔS(x), and S norm (x, θ i ) is the intensity value after normalization of the i-th angular projection data. T is a preset threshold used to filter points with low normalized intensity, thereby reducing the influence of background noise;

[0067] The specific process of compensation and calibration in step S65 is as follows:

[0068] S651: If Fill background data of size (Δ, θ) on the right side of the sinogram. If Fill background data of size (, θ) on the left side of the sinogram. The compensated sinogram The center is aligned with the projection rotation center. Among them, (Δ, θ) is the size of the background data to be filled. Δ represents the number of detectors to be filled, and θ represents all angles. is the sinogram after filling compensation, and the center has been aligned with the ideal rotation center, where is the compensated detector position;

[0069] S652: When the compensated is odd, Simultaneously cut off on both the left and right sides When the compensated is even, it will be adjusted according to The intensities of the boundary points on both sides are selected, and the stronger side is cut off The weaker side is cut off To achieve automatic calibration of the rotation center.

[0070] The beneficial effects of the present invention include:

[0071] The pulse - type terahertz CT data pre - processing method with adaptive recognition provided by the present invention collects multiple time - domain pulse signals with slight differences through a round - trip scanning mode at each fixed coordinate position, and performs averaging processing on the time - domain pulse signals of multiple repeated measurements at the same coordinate position; extracts the peak - to - peak value of the signal based on the peak - to - peak value extraction method with an adaptive threshold; detects the angular difference between each coordinate point and its adjacent points, selects the coordinate points that meet the interpolation conditions, and applies the linear interpolation method between the coordinate points that meet the interpolation conditions to generate uniformly distributed interpolation coordinates to achieve signal resampling; normalizes the peak - to - peak value of the collected signal and converts the signal from the peak - to - peak value to the absorption coefficient; performs flipping processing on even rows; realizes automatic calibration through a data - matching method and adjusts the rotation center of the sinogram. Through multiple samplings and averaging processing, the signal - to - noise ratio is improved, and the recognizability and stability of the signal are enhanced. The adaptive threshold method is used to accurately extract the peak - to - peak value of the signal to ensure a reliable characterization of the energy attenuation characteristics. Aiming at the non - uniformity in the sampling process, the present invention introduces a linear interpolation algorithm to unify the time and space scales and eliminate the acquisition error. In the signal characteristic analysis, based on the Lambert - Beer law, the absorption coefficient is calculated and combined with a threshold correction mechanism, effectively enhancing the processing robustness of low - signal - to - noise ratio data.

[0072] First, by performing multiple samplings and averaging processing on the signal, the signal - to - noise ratio is significantly improved, and the recognizability and stability of the signal are enhanced.

[0073] Second, the adaptive threshold method is used to accurately extract the peak - to - peak value of the signal to ensure a reliable characterization of the energy attenuation characteristics. Aiming at the non - uniformity in the sampling process, the present invention introduces a linear interpolation algorithm to unify the time and space scales and eliminate the acquisition error.

[0074] Third, in the signal characteristic analysis, based on the Lambert - Beer law, the absorption coefficient is calculated and combined with a threshold correction mechanism, effectively enhancing the processing robustness of low - signal - to - noise ratio data. In addition, to solve the problem of inconsistent data directions caused by round - trip scanning, the even - numbered rows of data are flipped to realize the normalization processing of the data format.

[0075] Finally, through edge detection, boundary point matching, and offset compensation, the rotation center of the sinogram is automatically calibrated, ensuring the accuracy and symmetry of tomographic image reconstruction, and providing a solid data foundation and processing framework for high - precision THz - CT applications. Brief Description of the Drawings

[0076] Figure 1 It is a flow schematic diagram of the pulse - type terahertz CT data pre - processing method with adaptive recognition of the present invention.

[0077] Figure 2 It is a schematic diagram of peak - to - peak value extraction of the present invention.

[0078] Figure 3 Schematic diagram of interpolating non-uniform data according to the present invention. Detailed implementation manners

[0079] The following will further elaborate on the present invention in conjunction with the Figures 1-3 accompanying drawings:

[0080] Embodiment 1

[0081] Refer to Figure 1 As shown, an adaptive recognition pulsed terahertz CT data preprocessing method includes the following steps:

[0082] S1: At each fixed (θ, y) coordinate position, multiple time-domain pulsed signals with slight differences are collected through a round-trip scanning mode, and the time-domain pulsed signals of multiple repeated measurements at the same coordinate position are averaged. The slight differences in the time-domain pulsed signals are mainly caused by unstable factors such as movement. Through the averaging process, the data volume can be effectively reduced, and at the same time, the signal-to-noise ratio can be significantly improved.

[0083] S2: Extract the peak-to-peak value of the signal based on the peak-to-peak extraction method with an adaptive threshold, which can ensure that the extracted peak-to-peak value can accurately reflect the energy attenuation characteristics of the signal;

[0084] S3: Detect the angular difference between each coordinate point and its adjacent points, select the coordinate points that meet the interpolation conditions, and apply the linear interpolation method between the coordinate points that meet the interpolation conditions to generate uniformly distributed interpolation coordinates, thereby realizing the resampling of the signal. Due to the influence of the round-trip movement of the sample during the scanning process, the collected coordinate data has velocity non-uniformity, which further leads to the inconsistency of the time and space scales of signal sampling. Therefore, an interpolation algorithm is used to process the collected coordinate data to make the time and space scales of signal sampling consistent.

[0085] S4: Normalize the peak-to-peak value of the collected signal, convert the signal from the peak-to-peak value to the absorption coefficient, calculate and correct the absorption coefficient of the signal, and quantitatively characterize the absorption characteristics of the sample;

[0086] S5: Flip the even rows, that is, the data corresponding to the reverse scan. Since the round-trip scanning mode is used for data collection, the sample moves back and forth along the y-axis between y max and y min during the scanning process. When the sample completes one scanning direction, it rotates by an angle for the next scan, and at this time, the direction of the collected data will be reversed, from y min to y max . To ensure the consistency of the data under different scanning directions, it is necessary to flip the even rows (the data corresponding to the reverse scan).

[0087] S6: Implement automatic calibration by means of data matching to adjust the rotation center of the sinogram.

[0088] Since terahertz CT is also affected by equipment and environmental noise as well as unstable factors during data acquisition, all of which may reduce the accuracy of the imaging results. Therefore, in order to overcome the limitations brought by these challenges and noise interference, the data preprocessing of the present invention is particularly important in the terahertz CT system. By optimizing the signal, correcting errors and suppressing noise, more reliable and accurate input data is provided for image reconstruction.

[0089] In the prior art, by collecting intermediate-frequency signals S(t) at different projection angles and azimuths and removing unnecessary background noise, background-removed intermediate-frequency signal data is obtained; the background-removed intermediate-frequency signal data is filtered to enhance the signal quality; using the phase gradient autofocus algorithm, point-by-point self-calibration is performed on the filtered intermediate-frequency signal data at different projection angles and azimuths to eliminate the phase error caused by equipment errors and environmental factors, and the self-calibrated intermediate-frequency signal S′(t) is obtained; the self-calibrated intermediate-frequency signal S′(t) of the frequency component is first subjected to Fourier transform and then peak extraction to obtain the projection data of the object This process has low processing efficiency, is prone to introducing human errors, and the accuracy and reliability of the processing results are insufficient.

[0090] In this embodiment, in the terahertz CT data preprocessing stage, by performing multiple samplings and averaging the signals, the signal-to-noise ratio is significantly improved, and the recognizability and stability of the signals are enhanced. At the same time, the adaptive threshold method is used to accurately extract the peak-to-peak value of the signal to ensure the reliable characterization of the energy attenuation characteristics. Aiming at the non-uniformity in the sampling process, the present invention introduces a linear interpolation algorithm to unify the time and space scales and eliminate the acquisition error. In the signal characteristic analysis, based on the Lambert-Beer law, the absorption coefficient is calculated and combined with the threshold correction mechanism to effectively enhance the processing robustness of low signal-to-noise ratio data. In addition, to solve the problem of inconsistent data directions caused by round-trip scanning, the even-numbered row data is flipped to realize the normalization processing of the data format. Finally, through edge detection, boundary point matching and offset compensation, the rotation center of the sinogram is automatically calibrated, ensuring the accuracy and symmetry of the tomographic image reconstruction. These optimization measures provide a solid data foundation and processing framework for high-precision THz-CT applications.

[0091] Embodiment 2

[0092] On the basis of Embodiment 1, the specific formula for averaging the time-domain pulse signals of multiple repeated measurements at the same coordinate position in step S1 is as follows:

[0093]

[0094] Among them, is the averaged signal, k is the total number of averaged signals, i is the single signal number, and S i is a single signal.

[0095] Through the above averaging process, noise interference can be effectively reduced and signal quality can be improved. Not only is the amount of data effectively reduced, but also the signal-to-noise ratio of the signal is significantly improved.

[0096] The specific process of extracting the peak-to-peak value of the signal based on the peak-to-peak extraction method with an adaptive threshold in step S2 is as follows:

[0097] S21: Extract each time-domain signal S i (t), and calculate the maximum and minimum values of each time-domain signal S i (t) to initially obtain the peak-to-peak value of the signal. The specific formulas for calculating the maximum and minimum values of the time-domain signal S i (t) are as follows:

[0098] max i = max(S i (t));

[0099] min i = min(S i (t));

[0100] Among them, max(·) is the maximum value function, and min(·) is the minimum value function;

[0101] S22: Introduce a vertical amplitude threshold to limit the peak amplitude, and combine a horizontal distance threshold to constrain the time interval between peaks;

[0102] S23: For signals exceeding the threshold limit, use the method of local extreme search to optimize the peak-to-peak extraction process, and optimize the calculation of the peak-to-peak value by finding local maximum and local minimum values within a window range;

[0103] S24: Define search windows Window max and Window min . The size of the window is determined by the horizontal distance threshold threshold x . Within the window Window max , find the local minimum local min , and within the window Window min , find the local maximum local max to obtain the final peak-to-peak value.

[0104] In this embodiment, the specific process of introducing the vertical amplitude threshold in step S22 to limit the peak amplitude and combining the horizontal distance threshold to constrain the time interval between peaks is as follows:

[0105] When the absolute value exceeds the vertical amplitude threshold, i.e. |max i |>threshold y , or |min i |>threshold y , then continue to calculate the peak-to-peak value. If it does not exceed, it is regarded as the background noise of the signal. The horizontal distance threshold constraint is threshold x ;

[0106] If the maximum value t max and minimum value t min The index difference is less than the horizontal distance threshold threshold x , that is |t max -t min | <threshold x , it is considered to be within the horizontal distance threshold x The lower maximum and minimum values come from a pulse signal, and the peak-to-peak value is calculated directly, Signal_FF=max i -min i .

[0107] In step S24, the search window is defined max and Window min The specific formula is as follows:

[0108] Window min =[t min -(threshold x -1), t min +(threshold x -1)];

[0109] Window min =[t min -(threshold x -1), t min +(threshold x -1)];

[0110] Finding local minimum min , and in the window Window min Find the local maximum max The specific formula for calculating the final peak-to-peak value is as follows:

[0111] Signal_FF = max(|max i | + |local min |, |min i | + |local max |);

[0112] Where, Window max is the maximum value of the search window on the abscissa time axis, Window min is the minimum value of the search window on the abscissa time axis, Signal_FF is the peak-to-peak value of the time-domain pulse signal, max i is the maximum value of the ordinate within the search window, min i is the maximum value of the ordinate within the search window, t max is the abscissa corresponding to max i , t min is the abscissa corresponding to min i .

[0113] See Figure 2 shown. To ensure that the extracted peak-to-peak value can accurately reflect the energy attenuation characteristics of the signal, this study adopted a peak-to-peak value extraction method based on an adaptive threshold. For each time-domain signal s i (t), first calculate its maximum and minimum values to initially obtain the peak-to-peak value of the signal. At the same time, introduce a vertical amplitude threshold to limit the peak amplitude, and combine a horizontal distance threshold to constrain the time interval between peaks. Vertical amplitude threshold constraint: If the absolute value exceeds the vertical amplitude threshold, that is, then continue to calculate the peak-to-peak value. If it does not exceed, it is regarded as the background noise of the signal. The horizontal distance threshold constraint is threshold x . If the index difference between the maximum value t max and the minimum value t min is less than the horizontal distance threshold. At this time, it is considered that the maximum and minimum values come from a pulse signal under the horizontal distance threshold threshold x , and then directly calculate the peak-to-peak value. For signals exceeding the threshold limit, further adopt a local extreme value search method to optimize the peak-to-peak value extraction process. At this time, the maximum value max i and the minimum value min i belong to different pulse signals respectively, and the peak-to-peak value calculation can be optimized by finding local maximum and local minimum values within a window range.

[0114] The local extreme value search not only optimizes the extraction of peaks, but also avoids the interference of signal noise on the peak-to-peak value calculation, ensuring the accurate reflection of the energy attenuation characteristics. In the signal, if the interval between the maximum value and the minimum value exceeds the horizontal distance threshold threshold x, they are considered different pulse signals. Depending on the requirements, pulse peaks at different time coordinates or the adaptive maximum peak value can be selected. This ensures reliable characterization of signal energy variations even under interference and noise conditions. This method enhances the robustness and accuracy of peak-to-peak value calculations, providing a reliable foundation for analyzing energy decay information in time-domain signals.

[0115] Example 3

[0116] On the basis of Example 1 or Example 2, in step S3, the specific process of applying the linear interpolation method between the coordinate points that meet the interpolation conditions to generate uniformly distributed interpolation coordinates is as follows:

[0117] S31: For each angle of the non-uniform coordinate x and the corresponding signal y, the maximum travel in the y direction is y max , first generate equally spaced horizontal coordinates x p , the horizontal axis is in the interval [-y max ,y max ] The number of points is controlled by the resolution r, that is,

[0118]

[0119] in, is the interval between the horizontal axes, and is the horizontal coordinate value that needs to be interpolated, k is the coordinate number corresponding to the coordinate that needs to be interpolated, and αx is the equal distance between the coordinates that need to be interpolated;

[0120] S32: Use linear interpolation to calculate the corresponding signal value y p , the specific formula is as follows:

[0121]

[0122] in, Is the difference function, through the input horizontal coordinate Calculate the corresponding signal value interpolation is the interpolation function, and linear is the linear interpolation function.

[0123] Due to the influence of the back-and-forth movement of the sample during the scanning process, the acquired coordinate data has velocity non-uniformity, which in turn causes inconsistencies in the time and space scales of signal sampling. To solve this problem, an interpolation algorithm is used to process the acquired coordinate data. By generating uniformly distributed interpolation coordinates, the original signal data can be converted into a signal under a uniform time scale, thereby eliminating the errors caused by scanning non-uniformity. Finally, the interpolated signal data not only has consistency in space but also ensures the uniformity of the time scale, providing a reliable data basis for subsequent image reconstruction and analysis.

[0124] The specific process of calculating and correcting the absorption coefficient of the signal in step S4 is as follows:

[0125] S41: Normalize the peak-to-peak value of the acquired signal to obtain the transmittance. The specific calculation formula is as follows:

[0126]

[0127] S42: Calculate the absorbance of the signal. The specific formula is as follows:

[0128]

[0129] Where A i is the absorbance of the i-th signal point, Si is the peak-to-peak value of the i-th signal, and s′ i is the normalized peak-to-peak value;

[0130] S43: Define the correction factor ∈. When A i is less than ∈, force its value to be not lower than ∈, that is, A i = max(A i , ∈).

[0131] Example 4

[0132] Based on Example 1 or Example 2 or Example 3, the specific process of flipping the data of the even rows in step S5, that is, the data corresponding to the reverse scan, is as follows:

[0133] By traversing the sinogram matrix, for each row of data starting from the second row of data, that is, the data with index 1, perform an element-by-element inversion operation so that the data direction of the reverse scan is consistent with the forward scan data of the previous row.

[0134] Since the data is acquired using a back-and-forth scanning mode, during the scanning process, the sample moves back and forth along the y-axis between y max and y min . When the sample completes one scanning direction, the rotation angle is adjusted for the next scan. At this time, the data acquisition direction will be reversed, from y min to ymax . To ensure the consistency of data under different scanning directions, it is necessary to perform a flipping process on the even rows (corresponding to the data of reverse scanning). The specific implementation method is to traverse the sinogram matrix, and for each row of data starting from the second row (index 1), perform an element-wise inversion operation, so that the data direction of reverse scanning is consistent with the forward scanning data of the previous row. This process effectively unifies the data format, eliminates the directional differences introduced by round-trip scanning, and provides a standardized data input for subsequent tomographic reconstruction and signal processing.

[0135] In step S6, automatic calibration is achieved through a data matching method, and the specific process of adjusting the rotation center of the sinogram is as follows:

[0136] S61: Let the data of the given sinogram be S(x, θ), where x represents the position of the detector and θ represents the projection angle; for each row of data, use the difference filter D = [1, -1] to calculate the difference of each row, and the specific formula is as follows:

[0137] ΔS(x) = |D * S(x, θ)|;

[0138] where * is the convolution operation, and ΔS(x) represents the change rate between two adjacent data;

[0139] S62: Normalize the intensity value of each row of data to extract the boundary points. The normalization of the i-th row of data can be expressed as:

[0140]

[0141] where ∈ is a specified constant, S(x, θ i ) is the projection data of the i-th angle, S min (x, θ i ) is the minimum value of the projection data of the i-th angle, S max (x, θ i ) is the maximum value of the projection data of the i-th angle;

[0142] S63: By calculating the change rate ΔS(x) and the normalized intensity S norm (x, θ i ), for each row of projection data, extract the points of the first and last changes in the data as the left boundary point x0 and the right boundary point x1;

[0143] S64: After extracting the left and right boundary points from each row, in order to minimize the deviation of the projection rotation center while retaining the position information of the measurement object, arrange the left boundary point x0 and the right boundary point x1 in ascending and descending order respectively, and take the mean of the first boundary points as the optimal boundary point for matching and

[0144] S65: Calculate the offset Δ between the left and right boundary points. By calculating the offset Δ, background data is filled on both sides of the sinogram for compensation, and the intensity of the specified boundary points is calibrated. The formula for calculating the offset Δ is as follows:

[0145]

[0146] In this embodiment, the specific formulas for extracting the first and last changed points in the data in step S63 as the left boundary point x and the right boundary point x1 are as follows:

[0147] x0 = arg min x (ΔS(x)>Q3(ΔS(x)) ∩ s norm (x, θ i )>T);

[0148] x1 = arg max x (ΔS(x)>Q3(ΔS(x)) ∩ S norm (x, θ i )>T);

[0149] Among them, argmin x returns the position where x is the smallest and satisfies the condition, argmax x returns the position where x is the largest and satisfies the condition, ΔS(x) is the change rate of two adjacent data, Q3 is the third quartile of ΔS(x), S norm (x, θ i ) is the normalized intensity value of the i-th angular projection data, and T is a preset threshold used to filter points with relatively low normalized intensity, thereby reducing the influence of background noise.

[0150] The specific process of compensation and calibration in step S65 is as follows:

[0151] S651: If fill background data of size (Δ, θ) on the right side of the sinogram, if fill background data of size (Δ, θ) on the left side of the sinogram. The compensated sinogram has its center aligned with the projection rotation center. Among them, (Δ, θ) is the size of the background data to be filled, Δ represents the number of detectors to be filled, and θ represents all angles. is the sinogram after filling and compensation, and its center has been aligned with the ideal rotation center, where is the compensated detector position;

[0152] S652: When the compensated is odd, Cut off both the left and right sides simultaneously When the compensated is an even number, it will be adjusted according to the intensities of the boundary points on both sides, and the stronger side will be selected to be cut off The weaker one is cut off Realize the automatic calibration of the rotation center. Through the above method, the automatic calibration of the rotation center is realized, ensuring the symmetry and matching of the data on both sides of the sinogram. This algorithm effectively solves the problem of reconstruction error caused by the offset of the rotation center, providing a reliable basis for high-precision CT image reconstruction.

[0153] In computer tomography (CT) data processing, to ensure the accuracy of image reconstruction, the calibration of the rotation center is crucial. When the rotation center of the collected sinogram data is not completely aligned with the geometric center of the object, it may lead to serious distortion of the reconstructed image. Therefore, it is necessary to achieve automatic calibration through data matching methods to adjust the rotation center of the sinogram.

[0154] In another implementation manner of this embodiment, the above automatic calibration method of the rotation center eliminates the offset of the rotation center caused by equipment errors, inconsistent scanning angles or other factors. To further improve the data quality and expand the adaptability to different scanning situations, the present invention covers a manual calibration method of the rotation center. This method provides flexibility based on the automatic calibration method and is applicable to scenarios where manual intervention is required to meet higher precision requirements in some special cases. Specifically, the manual calibration method observes the characteristics of the sinogram manually and makes targeted supplementary adjustments at the left and right ends of the data. This adjustment process is based on the visual analysis of the sinogram. After determining the area where data needs to be supplemented, the data is corrected by manually setting parameters to ensure the consistency between the rotation center and the geometric center of the object.

[0155] The present invention provides a preprocessing and correction method for pulsed terahertz computed tomography (THz-CT) data, characterized in that the method provides a systematic preprocessing and correction process. For the data collected by a terahertz time-domain spectroscopy computed tomography (THz-TDS-CT) system, given the path for data storage, through a series of automated processing steps, the projection data after noise suppression, signal calibration, and offset correction can be directly obtained. This provides reliable and accurate input data for subsequent image reconstruction. This process requires no manual intervention, greatly improving the processing efficiency and ensuring the accuracy and stability of image reconstruction. Through the extraction of the peak-to-peak value with an adaptive threshold, the signal peak value is accurately extracted, ensuring a reliable characterization of the signal energy attenuation characteristics. At the same time, the peak value calculation process is optimized, improving the robustness of signal processing and ensuring accuracy even under the influence of noise and unstable factors. The automatic calibration of the rotation center automatically calibrates the rotation center of the projection data through methods such as edge detection, boundary point matching, and offset compensation, thus avoiding reconstruction errors caused by the offset of the rotation center and improving the accuracy of image reconstruction.

[0156] In summary, the preprocessing method for pulsed terahertz CT data with adaptive recognition provided by the present invention collects multiple time-domain pulse signals with slight differences through a round-trip scanning mode at each fixed coordinate position, and performs averaging processing on the time-domain pulse signals of multiple repeated measurements at the same coordinate position; performs peak-to-peak extraction of the signal based on the peak-to-peak extraction method with an adaptive threshold; detects the angular difference between each coordinate point and its adjacent points, selects the coordinate points that meet the interpolation conditions, and applies a linear interpolation method between the coordinate points that meet the interpolation conditions for processing to generate uniformly distributed interpolation coordinates to achieve signal resampling; normalizes the peak-to-peak value of the collected signal and converts the signal from the peak-to-peak value to the absorption coefficient; performs flipping processing on the even rows; realizes automatic calibration through data matching methods to adjust the rotation center of the sinogram. Through multiple samplings and averaging processing, the signal-to-noise ratio is improved, enhancing the recognizability and stability of the signal. The adaptive threshold method is used to accurately extract the peak-to-peak value of the signal, ensuring a reliable characterization of the energy attenuation characteristics. Aiming at the non-uniformity in the sampling process, the present invention introduces a linear interpolation algorithm to unify the time and space scales and eliminate the acquisition error. In the signal characteristic analysis, based on the Lambert-Beer law, the absorption coefficient is calculated and combined with a threshold correction mechanism, effectively enhancing the processing robustness of low signal-to-noise ratio data.

[0157] By sampling multiple times and averaging the signals, the signal-to-noise ratio is significantly improved, and the recognizability and stability of the signals are enhanced. The adaptive threshold method is used to accurately extract the peak-to-peak value of the signals to ensure a reliable characterization of the energy attenuation characteristics. Aiming at the non-uniformity in the sampling process, the present invention introduces a linear interpolation algorithm to unify the time and space scales and eliminate the acquisition error. In the analysis of signal characteristics, based on the Lambert-Beer law, the absorption coefficient is calculated and combined with a threshold correction mechanism, effectively enhancing the processing robustness of low signal-to-noise ratio data. In addition, to solve the problem of inconsistent data directions caused by round-trip scanning, the even-row data is flipped to normalize the data format. Through edge detection, boundary point matching, and offset compensation, the rotation center of the sinogram is automatically calibrated, ensuring the accuracy and symmetry of tomographic image reconstruction, providing a solid data foundation and processing framework for high-precision THz-CT applications.

Claims

1. A pulse-type terahertz CT data preprocessing method for adaptive recognition, characterized in that It includes the following steps: S1: At each fixed coordinate position, collect multiple time-domain pulse signals with slight differences through a round-trip scanning mode, and perform averaging processing on the time-domain pulse signals obtained from multiple repeated measurements at the same coordinate position; S2: Extract the peak-to-peak value of the signal based on the peak-to-peak extraction method with an adaptive threshold; S3: Detect the angular difference between each coordinate point and its adjacent points, select the coordinate points that meet the interpolation conditions, and apply the linear interpolation method between the coordinate points that meet the interpolation conditions to generate uniformly distributed interpolation coordinates, realizing resampling of the signal; S4: Normalize the peak-to-peak value of the collected signal, convert the signal from the peak-to-peak value to the absorption coefficient, and calculate and correct the absorption coefficient of the signal; S5: For even rows, that is, the data corresponding to the reverse scan, perform flipping processing; S6: Achieve automatic calibration through the method of data matching, and adjust the rotation center of the sinogram.

2. The pulse-type terahertz CT data preprocessing method for adaptive recognition according to claim 1, characterized in that The specific formula for performing averaging processing on the time-domain pulse signals obtained from multiple repeated measurements at the same coordinate position in step S1 is as follows: Among them, is the averaged signal, k is the total number of averaged signals, i is the single signal number, and S i is a single signal.

3. An adaptive recognition-based pulse terahertz CT data preprocessing method according to claim 1, characterized in that, The specific process of extracting the peak-to-peak value of the signal based on the peak-to-peak extraction method with an adaptive threshold in step S2 is as follows: S21: Extract each time-domain signal S i (t), and calculate the maximum and minimum values of each time-domain signal S i (t) to initially obtain the peak-to-peak value of the signal, and calculate the specific formulas for the maximum and minimum values of the time-domain signal S i (t) are as follows: max i = max(S i (t)); min i = min(S i (t)); Among them, max(·) is the function for taking the maximum value, and min(·) is the function for taking the minimum value; S22: Introduce a vertical amplitude threshold to limit the peak amplitude, and combine it with a horizontal distance threshold to constrain the time interval between peaks; S23: For the signals that exceed the threshold limit, use the method of local extreme value search to optimize the process of peak-to-peak value extraction, and optimize the calculation of the peak-to-peak value by finding local maxima and local minima within a window range; S24: Define the search window Window max and Window min , the size of the window is determined by the horizontal distance threshold threshold x . Within the window Window max , find the local minimum local min , and find the local maximum local min within the window Window max to obtain the final peak-to-peak value.

4. An adaptive recognition-based pulsed terahertz CT data preprocessing method according to claim 3, characterized in that, The specific process of introducing a vertical amplitude threshold to limit the peak amplitude and combining it with a horizontal distance threshold to constrain the time interval between peaks in step S22 is as follows: When the absolute value exceeds the vertical amplitude threshold, i.e., |max i | > threshold y , or |min i | > threshold y , then continue to calculate the peak-to-peak value. If it does not exceed, it is regarded as the background noise of the signal; the horizontal distance threshold is constrained to threshold x ; If the maximum value t max and the minimum value t min have an index difference less than the horizontal distance threshold threshold x , that is, |t max -t min |<threshold x , then it is determined that the maximum and minimum values come from a pulse signal under the horizontal distance threshold threshold x , and the peak-to-peak value is directly calculated, Signal_FF = max i -min i .

5. The pulse-type terahertz CT data preprocessing method with adaptive recognition according to claim 3, characterized in that Define the search window Window in step S24 max and Window min The specific formula is as follows: Window max = [t max -(threshold x - 1), t max +(threshold x - 1)]; Window min = [t min -(threshold x - 1), t min +(threshold x - 1)]; Search for local minimum min and search for local maximum within the window Window min The specific formula for calculating the final peak-to-peak value after that is as follows: max ​ Signal_FF = max(|max i | + |local min |, |min i | + |local max |); Among them, Window max is the maximum value of the search window on the horizontal time axis, Window min is the minimum value of the search window on the horizontal time axis, Signal_FF is the peak-to-peak value of the time-domain pulse signal, max i is the maximum value of the vertical coordinate within the search window, min i is the maximum value of the vertical coordinate within the search window, t max is for max i corresponding abscissa, t min is for min i corresponding abscissa.

6. An adaptive recognition-based pulsed terahertz CT data preprocessing method according to claim 1, characterized in that The specific process of applying the linear interpolation method between the coordinate points that meet the interpolation conditions to generate uniformly distributed interpolation coordinates in step S3 is as follows: S31: For the uneven coordinates x and the corresponding signals y at each angle, the maximum travel in the y direction is y max , first generate equally spaced abscissas x p , where the abscissas are equally spaced values in the interval [-y max , y max , and the number of points is controlled by the resolution r, i.e., Among them, is the interval between abscissas, and is the abscissa value to be interpolated, k is the number of coordinates corresponding to the coordinate to be interpolated, and Δx is the equal interval distance of the coordinate to be interpolated; S32: Use linear interpolation to calculate the corresponding signal value The specific formula is as follows: Among them, is the difference function, and the corresponding signal value is calculated through the input abscissa interpolation is the interpolation function, and linear is the linear interpolation function.​ 7. An adaptive recognition-based pulse terahertz CT data preprocessing method according to claim 1, characterized in that The specific process of calculating and correcting the absorption coefficient of the signal in step S4 is as follows: S41: Normalize the peak-to-peak value of the collected signal to obtain the transmittance, and the specific calculation formula is as follows: S42: Calculate the absorbance of the signal, and the specific formula is as follows: Among them, A i is the absorbance of the i-th signal point, Si is the peak-to-peak value of the i-th signal, and s i ′ is the peak-to-peak value after normalization; S43: Define a correction factor ∈. When A i is less than ∈, force its value to be not lower than ∈, that is, A i = max(A i , ∈).

8. An adaptive recognition-based pulsed terahertz CT data preprocessing method according to claim 1, characterized in that, The specific process of performing flipping processing on even rows, that is, the data corresponding to the reverse scan, in step S5 is as follows: By traversing the sinogram matrix, for each row of data starting from the second row of data, that is, the data with an index of 1, perform an element-by-element inversion operation, so that the data direction of the reverse scan is consistent with the forward scan data of the previous row.

9. An adaptive recognition-based pulsed terahertz CT data preprocessing method according to claim 1, characterized in that The specific process of achieving automatic calibration through the method of data matching and adjusting the rotation center of the sinogram in step S6 is as follows: S61: Let the data of the given sinogram be S(x,θ), where x represents the position of the detector and θ represents the projection angle; for each row of data, use the difference filter D = [1, -1] to calculate the difference of each row, and the specific formula is as follows: ΔS(x) = |D * S(x,θ)|; Among them, * is the convolution operation, and ΔS(x) represents the change rate between two adjacent data; S62: Normalize the intensity values of each row of data to extract boundary points. The normalization of the i-th row of data can be expressed as: where the correction factor ∈ is a specified constant, S(x,θ i ) is the projection data at the i-th angle, S min (x,θ i ) is the minimum value of the projection data at the i-th angle, S max (x,θ i ) is the maximum value of the projection data at the i-th angle; S63: By calculating the change rate ΔS(x) and the normalized intensity S norm (x, θ i ), for each row of projection data, extract the points where the data first and last change as the left boundary point x0 and the right boundary point x1; S64: After extracting the left and right boundary points from each line, to minimize the deviation of the projection rotation center while retaining the position information of the measurement object, the left boundary point x0 and the right boundary point x1 are sorted in ascending and descending order respectively, and the mean of the first boundary points is taken as the optimal boundary point for matching and S65: Calculate the offset Δ between the left and right boundary points. By calculating the offset Δ, background data is filled on both sides of the sinogram for compensation, and the intensity of the specified boundary points is calibrated. The formula for calculating the offset Δ is as follows:

10. A pulse terahertz CT data preprocessing method for adaptive recognition according to claim 9, characterized in that, The specific formulas for extracting the points of the first and last changes in the data as the left boundary point x0 and the right boundary point x1 in step S63 are as follows: x0 = arg min x (ΔS(x) > Q3(ΔS(x)) ∩ S norm (x, θ i )) > T); x1 = arg max x (ΔS(s)>Q3(ΔS(x)) ∩ S norm (x, θ i ))>T); where argmin x returns the position where x is the smallest and satisfies the condition, and argmax x returns the position where x is the largest and satisfies the condition. ΔS(x) is the change rate between two adjacent data, Q3 is the third quartile of ΔS(x), and S norm (x, θ i ) is the intensity value after normalization of the i-th angular projection data. T is a preset threshold used to filter points with low normalized intensity, thereby reducing the influence of background noise; The specific processes of compensation and calibration in step S65 are as follows: S651: If background data of size (Δ, θ) is filled on the right side of the sinogram, if background data of size (Δ, θ) is filled on the left side of the sinogram, the compensated sinogram has its center aligned with the projection rotation center, where (Δ, θ) is the size of the background data to be filled, Δ represents the number of detectors to be filled, and θ represents all angles, is the sinogram after filling compensation, and its center has been aligned with the ideal rotation center, where is the compensated detector position; S652: When the compensated is odd, simultaneously cut off from both the left and right sides When the compensated is even, it will be adjusted according to the intensities of the boundary points on both sides, and the stronger side will be selected for cutting off the weaker one is cut off to achieve automatic calibration of the rotation center.

Images