CN114283220A - Image reconstruction and system error correction method for cone CT - Google Patents

Abstract

The invention discloses a method for image reconstruction and system error correction of cone CT. The method comprises the steps of firstly carrying out CT scanning on a sample to obtain original data and simultaneously carrying out bright/dark current correction, quickly obtaining a reconstructed CT original image by utilizing the cooperative operation of a CUDA and a multi-thread CPU according to the corrected data, reducing the total time while reducing the requirements on a memory and a video memory by segmented reconstruction streaming transmission, calculating a ray source optical axis deviation value, a detector center deviation value, an object distance and a detector distance of an imaging system according to the two-bit gray level entropy of a tomographic image at a specific position, and correcting the parameters and obtaining a final CT reconstructed image during back projection reconstruction. Compared with the prior art, the method has the advantages that the speed and the quality are improved, meanwhile, the requirements on the operation space and the precision of each component of the CT imaging system are reduced, and meanwhile, the calculation of the space parameters is completed under the condition that a calibration object is not used.

Description

Image reconstruction and system error correction method for cone CT

Technical Field

The invention relates to the field of CT imaging, in particular to a method for image reconstruction and system error correction of cone CT.

Background

Cone CT reconstruction is typically very time consuming, requires relatively many computation spaces and computation performance, and requires relatively high precision algorithms for each part of the imaging system. For example, in a 1500 × 2000 pixel detector, 720 projection views are acquired by rotating one circle, and when the detector distance is 10000 times the detector pixel side length, the original data itself occupies a memory larger than 4GB, the CUDA acceleration needs to occupy a video memory larger than 8GB, and the reconstruction result occupies about 10 to 30GB according to the difference of the reconstruction field of view and the resolution. Meanwhile, in the aspect of system errors, under the condition that the center of the rotating shaft and the center of the detector have only 1-degree deviation, through the amplification of the conical CT light path, deviation of more than 15 pixel points can be caused on each projection drawing finally, and the final result is greatly influenced.

The existing solution to the operational problem is to add a memory and a video memory, and to allocate a memory above 100GB to process all projection images at one time. After the hardware conditions are determined, the upper limit of the reconstruction precision of the whole system is determined simultaneously by the schemes, and the reconstruction exceeding the design precision of the schemes cannot be carried out. In addition, the steps of loading the projection diagram, transmitting source data to the GPU, performing fault diagram reconstruction operation on the GPU, reading a result from the GPU and the like are relatively independent, which causes the GPU to be in an idle state when data is transmitted, and causes waste of computing power when the GPU is in an idle state when data is transmitted.

In the aspect of solving the system error, the existing scheme usually uses hardware means, including measuring the distance of the rotating table by a grating ruler, aligning and adjusting the position of the optical axis in advance by an optical sighting device, additionally installing a marble base to reduce the deformation caused by the temperature, and the like. These approaches are costly and can significantly increase the cost of the CT system. Meanwhile, when the system error is increased due to some reason, only the quality reduction of the reconstruction result can be observed, the error reason can not be detected through software, and only hardware inspection can be carried out.

Therefore, there is a need to provide a new CT reconstruction method to solve the problems of the prior art.

Disclosure of Invention

Technical problem to be solved

The invention aims to solve the technical problem that: how to accelerate the CT reconstruction speed, reduce the performance requirement while maintaining the reconstruction image precision, and calculate the system hardware error through a software algorithm.

(II) technical scheme

In order to solve the above problems, the present invention provides the following technical solutions, which provide a method for image reconstruction and system error correction for cone CT, and the method is as follows.

A method for image reconstruction and systematic error correction for cone CT, comprising the steps of:

a, carrying out CT scanning on a sample to obtain original data;

b, reconstructing a CT original image by utilizing the flow line type cooperative operation of the CUDA and the multi-thread CPU according to the corrected data;

step C, calculating hardware error parameters of the imaging system according to the two-position gray entropy of the tomogram at the middle position and the tomograms at the higher and lower positions, wherein the hardware error parameters specifically comprise a rotation center deviation value, a detector center deviation value, an object distance and a detector distance of the imaging system;

and D, correcting hardware error parameters during back projection reconstruction to obtain a final CT reconstructed image.

Preferably, the pipelined cooperation used in step B specifically comprises the steps of:

step B1, the CPU recombines the original CT projection drawing into a sine-wave drawing set by taking the size of the system memory as the limit;

step B2, transmitting the recombined sine pattern set to a GPU video memory by taking the video memory size as a limit;

b3, selecting according to actual conditions, calculating a new sine map set, or adjusting a reconstruction task responsible by part of the GPU to a CPU for processing;

step B4, the GPU reconstructs a fault map based on the sinogram in the display memory and informs the CPU which sinograms are not needed any more;

a step B5 of transmitting the tomogram result to the analysis computer through the network;

the steps B1, B2 and B3 are handled by the CPU, and the steps B4 and B5 are handled by the GPU, which are relatively independent.

Preferably, in the step C, a hardware error is calculated by a software means based on the tomogram, and the specific steps are as follows:

step C1, reconstructing the tomogram at the midpoint of the target area separately, and reconstructing all the median tomograms with the offset value of the detector in the x direction ranging from-10 mm to 10mm at intervals of 0.1 mm; selecting the clearest picture from the two-dimensional entropy as a standard, and recording a corresponding detector x-direction deviation value;

step C2, reconstructing the median tomogram independently, substituting the offset value in the x direction, and obtaining the offset value in the y direction of the detector by the same method standard;

step C3, reconstructing a higher fault and a lower fault by using FBP algorithm, calculating the x deviation of the two layers by using the method in step C1, wherein the difference value of the two deviations is the deviation value in the x direction of the rotation center, namely the higher the x-axis deviation of the position caused by the rotation axis not strictly vertically upwards is; the x deviation value outside the two layers is obtained by linear interpolation;

step C4, a rough detector-ray source distance and a rotation axis-ray source distance are given, an FDK algorithm is used for reconstructing a higher fault and a lower fault, 1mm is taken as an interval, all fault graphs with a rotation axis-ray source distance deviation value ranging from-100 mm to 100mm are reconstructed, the clearest picture is selected by taking the lowest two-dimensional entropy as a standard, and the corresponding rotation axis-ray source distance is recorded;

step C5, carrying out global stepping optimization on the basis of the parameters, taking 0.1mm as a step length, adding small changes to each system hardware parameter including a rotation center deviation value, a detector center deviation value, an object distance and a detector distance, and investigating whether the picture is clear or fuzzy, wherein a parameter combination when the picture is clearest is taken as an optimal parameter combination;

and step C6, obtaining final parameters by adopting the steps, and reconstructing all the tomograms.

Preferably, in the error calculation process, the tomogram is segmented by using an OTSU algorithm, and two-dimensional entropy is calculated for the segmented image to evaluate the definition of the reconstructed image.

(III) advantageous effects

Compared with the prior art, the method has the advantages of obvious and positive technical effects and is mainly divided into two aspects. Firstly, on the basis of the existing technology of accelerating reconstruction by using the GPU, the invention improves the algorithm flow, so that the CPU and the GPU can better cooperate with each other, the waste of computing power caused by idle of a certain party is avoided, and compared with the scheme of reconstructing by using only the GPU, the reconstruction speed is optimized in time. And different from the route of integrally transmitting, integrally calculating and integrally obtaining the result in other schemes, the technical scheme of the invention adopts a pipeline form, reduces the consumption in the calculation process in space and can reconstruct the reconstruction domain exceeding the maximum memory space of the system. Secondly, different from the two methods for obtaining hardware parameters such as ray source optical axis offset, detector center offset, object distance, detector distance and the like in the prior art, namely, accurate hardware equipment is needed to reduce system errors and provide hardware parameters, or a specific calibration object is needed to perform hardware calibration once before actual shooting. The tolerance of the system to hardware errors can be improved, the manufacturing cost and the maintenance cost are reduced, and the settings such as the integral amplification factor of the system can be set by a user without extra hardware calibration.

Drawings

FIG. 1 is a flowchart of a CT reconstruction method according to the present invention

FIG. 2 is a schematic diagram of a CT image acquisition system

FIG. 3 is a CPU/GPU cooperative reconstruction process used in the present invention

FIG. 4 is a schematic diagram of the effect of reconstruction of a fault using different total angle numbers

FIG. 5 is a schematic diagram of the light path of a cone CT in three different situations

FIG. 6 is a schematic diagram of using OTSU algorithm to segment and measure the detector center offset by the result

Detailed Description

The invention is further explained below with reference to the drawings and the embodiments.

The method for image reconstruction and system error correction for cone CT of the present invention is shown in FIG. 1 and comprises the following steps.

Step 1, carrying out CT scanning on a sample to obtain original data.

The raw data is the signals output by the detectors in the CT device indicating the respective intensities of the radiation they receive, and if the CT device shown in fig. 2 is used, the raw projection data is cone beam projection data, and the signals are then represented as a two-dimensional projection of the sample. The total number of sheets is related to the turntable in fig. 2 and the currently set refresh rate of the detector.

And 2, reconstructing a CT original image by utilizing a flow line type cooperative operation of a CUDA and a multi-thread CPU according to the corrected data.

FIG. 3 is a schematic diagram of this step, which can more intuitively see the work that CPU and GPU respectively play, the CPU is responsible for reconstructing the sine map set and the GPU is responsible for reconstructing the sine map set into a fault map. The speed of preparing the sine atlas by the CPU is faster than the speed of reconstructing the fault by using the atlas by the GPU, and the size of the video memory is generally smaller than that of the memory, so that the video memory space is full for a period of time, the data in the video memory space is not processed, but the CPU is in an idle state. At the moment, the algorithm can evaluate the time of completing reconstruction by the GPU and the time of reconstructing a fault by using a CPU multithread, and bears a part of reconstruction tasks to accelerate the speed when the CPU is idle.

And 4, calculating system hardware parameters such as a rotation center deviation value, a detector center deviation value, an object distance, a detector distance and the like of the imaging system according to the two-bit gray entropy of the tomograms at specific positions such as the tomogram at the middle position, the tomograms at the higher and lower positions and the like.

The reconstruction algorithm used in the present invention is based on a filtered back-projection algorithm. Filtering is used to improve the reconstruction quality and backprojection is actually based on ray paths, summing the actual values of all angles of rays passing through a certain voxel point to obtain the density of the point. The more angles are used, the better the reconstructed image effect is. Figure 4 shows the effect of the tomogram when reconstructed using 4, 8, 15, 25, 50 and 1440 angles. As can be seen from the figure, an important step in reconstructing the map is to determine which pixel point in the projection map corresponds to the ray passing through a certain point P (x, y) in the back projection process of a certain angle.

First, consider the tomographic reconstruction at the midpoint elevation, i.e., on a horizontal plane including the optical axis. At the level of the midpoint, the rays do not travel up and down, but hit a line in the midpoint of the detector. The pixels at different positions of the row represent the attenuation rates of the rays at different positions.

Fig. 5 shows the imaging diagram at the level of the midpoint when the rotation table is rotated to a certain angle α. The left image is the ideal situation, i.e. the optical axis of the radiation source, the central axis of the rotation stage and the central point of the detector are all aligned exactly. In the figure, S is a ray source, O is a central point of a rotating shaft, and is also a coordinate origin and a line segment when the plane carries out fault reconstruction

And C is the central point of the detector, the optical axis of the ray from the ray source passes through the central point O of the rotating shaft and then intersects with the detector on C, and the point is also the central point of the detector.

The intensity corresponding to the point P with the coordinate (x, y) at the angle is the detector at P₀The signal value of (2) can be calculated from the trigonometric relation

The length of (A) is as follows:

wherein D_dDistance of the source to the detector, i.e. line segment in the figure

Length of (d). D_sDistance of the radiation source from the rotation center of the turntable, i.e. line segment in the figure

Length of (d). x and y are the coordinates of the voxel point being reconstructed at the current height, with the origin of coordinates at the center of rotation.

The density value of the voxel P (x, y) at this time is:

wherein Q_αWhen the rotation angle is alpha, the projection function after weighted and convolution filtering is carried out,

represents the abscissa position of the target pixel point in the projection drawing, and 0.5W is the middle position of the detector where the point C is located. 0.5H represents the height at which the target voxel is at the very center at this time.

The image in fig. 5 is an imaging schematic diagram in an actual situation, where the optical axis of the radiation source, the rotation axis of the rotation stage and the center point of the detector are not exactly aligned. Similarly, S in the figure is a ray source, O is a center point of a rotating shaft, and is also an origin of coordinates and a line segment when the plane is subjected to tomographic reconstruction

Representing a flat panel detector. At this time, a ray from the radiation source passing through the rotation axis center point O intersects the detector at a point C. And the optical axis of the ray source intersects with the detector at a point D₀The actual center point of the detector is located at point D.

For convenient operation, the ray passing through the center of the rotating shaft, namely the ray of the coordinate origin O in the process of reconstructing the tomogram

For reference, the optical axis of the ray source is opposite

Is offset from

Is recorded as the optical axis deviation O of the ray source_LOffset of center point D of detector with respect to C

Is recorded as the center offset O of the detector_D. Similarly, the intensity corresponding to the point P with coordinates (x, y) at this angle is the detector at P₀The signal value of (2) can be calculated from the trigonometric relation

The length of (A) is as follows:

in this formula, when O_LAnd O_DWhen the value is 0, the value is obtained in the formula (1)

The result of the length. Likewise, the density value of the voxel P (x, y) at this time is:

in the two cases, only the tomographic reconstruction at the middle height is considered, and the projections of all points in the same tomographic plane in the flat panel detector are located at the same height, namely in the same row. Now consider the reconstruction in other planes, and in general, points at different positions in one slice plane will project at different heights on the flat panel detector, since some are closer to the source and some are further from the source. The right diagram of fig. 5 shows the light path diagram at this time.

For a point P (x, y, h) with coordinates (x, y) on the slice plane with horizontal height h, the corresponding ray at angle α also hits the flat panel detector on the abscissa shown by equation (3), while the ordinate position is:

likewise, the density value of the voxel P (x, y) at this time is:

from the above formula, the final result is subject to the optical axis offset O of the radiation source_LOffset of the center of the detector O_DDistance D from source to detector_dAnd the distance D from the ray source to the rotation center of the rotating table_sThe influence of (c). Generally, these parameters need to be obtained from hardware devices or kept stable after being fixed, and the technical solution of the present invention provides a means for calculating these parameters by software algorithms.

Offset of optical axis of ray source O_LCan be rewritten as:

O_L＝D_d tan β (7)

wherein beta is the angle of the optical axis of the ray source offset relative to the center of rotation when D_dThe offset O of the optical axis of the ray source is increased when the flat panel detector is placed at a longer distance_LAnd will increase with it, then there are:

the formula is composed of two parts, one is a detectorCenter offset O_DOne is the distance D from the source to the detector_dA proportional value. This indicates that O_DThe other three quantities are not mutually influenced, and the result is consistent with the common knowledge. Can first pass the evaluation of different O_DThe quality of the reconstructed image obtained by value taking is inversely obtained to obtain an optimal O_DThe value is obtained.

FIG. 6 is a graph of the difference O_DThe fault reconstruction result obtained under the value taking is obtained, the first line is the actual reconstruction result, the second line is the result after OTSU segmentation, the graph with the least bright spots after segmentation is selected, and the graph with the least and most clear diffuse spots in all the graphs is obtained, so that the detector center offset O of the system used in the example can be obtained_DIs-12.

To adjust O completely_DThen, the second part is reacted with D_dIs proportional, meaning that D is here_dThe quality of the picture is not affected, but the size of the picture break layer part is affected. With other parameters fixed, the farther away the detector is, the larger the image will naturally be. Thus, a fixed D can be temporarily set at will_dValues to first determine the size of the remaining two parameters. Since β essentially belongs to the systematic error due to misalignment of the optical path, and D_sIs a parameter set by the user, so that β is usually relatively small in actual processing. D_sThe disturbance caused is usually larger than beta, so O is first determined and used here_DThe same method firstly obtains a D_sAnd then determines the value of β.

D_dDoes not affect the definition of the tomogram at the middle position, needs to select a proper height h, and measures D by using the tomogram at the height h_dThe value of (c). A height is generally chosen that is sufficiently far from the midpoint to encompass a sample of sufficiently large cross-sectional area. As can be seen from the equations (5) and (6), D_dThe value of (a) will greatly influence the pixel positions corresponding to the calculated projection map of the voxel points at different positions. The same method can be conveniently used for the measurement.

And 5, correcting the parameters during back projection reconstruction and obtaining a final CT reconstruction image.

The specific flow is shown in fig. 3, but is slightly different from the previous one. In the above steps, the slice reconstruction is only a planar slice reconstruction, corresponding sinograms can be loaded according to actual needs, and the last step is reconstruction of all slices, namely, all sinograms are loaded by taking video storage as a limit and then reconstruction is carried out together.

The specific examples described in the application are only illustrative of the spirit of the invention. Various modifications, additions and substitutions of types may be made by those skilled in the art without departing from the spirit of the invention or exceeding the scope of the invention as defined in the accompanying claims.

Claims (4)

1. A method for image reconstruction and systematic error correction for cone CT, comprising the steps of:

a, carrying out CT scanning on a sample to obtain original data;

b, reconstructing a CT original image by utilizing the flow line type cooperative operation of the CUDA and the multi-thread CPU according to the corrected data;

step C, calculating hardware error parameters of the imaging system according to the two-position gray entropy of the tomogram at the middle position and the tomograms at the higher and lower positions, wherein the hardware error parameters specifically comprise a rotation center deviation value, a detector center deviation value, an object distance and a detector distance of the imaging system;

and D, correcting hardware error parameters during back projection reconstruction to obtain a final CT reconstructed image.

2. The method for image reconstruction and systematic error correction for cone CT according to claim 1, wherein the pipelined cooperation used in step B comprises the following steps:

step B1, the CPU recombines the original CT projection drawing into a sine-wave drawing set by taking the size of the system memory as the limit;

step B2, transmitting the recombined sine pattern set to a GPU video memory by taking the video memory size as a limit;

b3, selecting according to actual conditions, calculating a new sine map set, or adjusting a reconstruction task responsible by part of the GPU to a CPU for processing;

step B4, the GPU reconstructs a fault map based on the sinogram in the display memory and informs the CPU which sinograms are not needed any more;

a step B5 of transmitting the tomogram result to the analysis computer through the network;

the steps B1, B2 and B3 are handled by the CPU, and the steps B4 and B5 are handled by the GPU, which are relatively independent.

3. The method for image reconstruction and systematic error correction for cone CT as claimed in claim 1, wherein in step C the hardware error is calculated by software means based on the tomogram, specifically comprising the steps of:

step C1, reconstructing the tomogram at the midpoint of the target area separately, and reconstructing all the median tomograms with the offset value of the detector in the x direction ranging from-10 mm to 10mm at intervals of 0.1 mm; selecting the clearest picture from the two-dimensional entropy as a standard, and recording a corresponding detector x-direction deviation value;

step C2, reconstructing the median tomogram independently, substituting the offset value in the x direction, and obtaining the offset value in the y direction of the detector by the same method standard;

step C3, reconstructing a higher fault and a lower fault by using FBP algorithm, calculating the x deviation of the two layers by using the method in step C1, wherein the difference value of the two deviations is the deviation value in the x direction of the rotation center, namely the higher the x-axis deviation of the position caused by the rotation axis not strictly vertically upwards is; the x deviation value outside the two layers is obtained by linear interpolation;

step C4, a rough detector-ray source distance and a rotation axis-ray source distance are given, an FDK algorithm is used for reconstructing a higher fault and a lower fault, 1mm is taken as an interval, all fault graphs with a rotation axis-ray source distance deviation value ranging from-100 mm to 100mm are reconstructed, the clearest picture is selected by taking the lowest two-dimensional entropy as a standard, and the corresponding rotation axis-ray source distance is recorded;

step C5, carrying out global stepping optimization on the basis of the parameters, taking 0.1mm as a step length, adding small changes to each system hardware parameter including a rotation center deviation value, a detector center deviation value, an object distance and a detector distance, and investigating whether the picture is clear or fuzzy, wherein a parameter combination when the picture is clearest is taken as an optimal parameter combination;

and step C6, obtaining final parameters by adopting the steps, and reconstructing all the tomograms.

4. The method for image reconstruction and systematic error correction for cone CT according to claim 3, characterized in that during the error calculation, the tomography image is segmented using OTSU algorithm, and the two-dimensional entropy is calculated for the segmented image to evaluate the sharpness of the reconstructed image.

Images