RU2816842C1 - System of automatic search and correction of object rotation axis - Google Patents

Abstract

FIELD: X-ray tomography.

SUBSTANCE: for automatic search and correction of the axis of rotation of an object when performing X-ray tomography. The essence of the invention is that a system for automatically searching and correcting the axis of rotation of an object when performing X-ray tomography contains an X-ray source, a device for rotating and positioning the object, a flat panel detector and a computing device, wherein the flat panel detector is configured to generate projection data of the object, and the computing device is configured to: a) obtain projection data, b) calculate the average projection image, c) preprocess the average projection image using linear normalization and scaling, d) two-stage search for the values of the angle of inclination and shift of the rotation axis, e) clarify the position of the object's rotation axis.

EFFECT: providing the opportunity to improve the quality of reconstruction of the internal structure of the object under study.

1 cl, 13 dwg, 1 tbl

Description

The invention relates to the field of X-ray engineering and can be used both in medicine and in the field of technical non-destructive testing.

A projection method for positioning the center of rotation of a computed tomography system is known from the prior art, disclosed in patent document CN 102539460 A, published on 07/04/2012. The invention discloses a projection method for positioning the center of rotation of a computed tomography system. The method includes the steps of: obtaining 360 degree projection data to obtain a projection sine wave diagram; setting a proper threshold, segmenting the projection sinusoidal diagram, overlaying and averaging the projection data of all projection angles according to the corresponding pixels, i.e. overlaying all lines of the projection sinusoidal diagram and averaging to obtain a one-dimensional signal y=f(s), using the quadratic polynomial y=as2+bs+c and least squares fitting y=f(s); and calculating the value of the horizontal coordinates s=b/(-2a) of the symmetrical center of the parabola y=as2+bs+c, where the values of the horizontal coordinates are the coordinates of the center of rotation of the projection.

The disadvantage of the prototype is the impossibility of using it for the cone model of X-ray propagation, which is currently the most common model.

The objective of the invention is to eliminate the shortcomings of the prototype.

The technical result of the invention is to improve the quality of reconstruction of the internal structure of the object under study.

The technical result is achieved due to the fact that the system for automatically searching and correcting the axis of rotation of an object when performing X-ray tomography contains an X-ray source, a device for rotating and positioning the object, a flat panel detector and a computing device, characterized in that the flat panel detector is configured to generate projection data of the object , and the computing device is configured to:

a) obtaining projection data,

b) calculating the average projection image,

c) preprocessing the average projection image using linear normalization and scaling,

d) two-stage search for the values of the angle of inclination and shift of the axis of rotation,

e) clarifying the position of the object’s rotation axis.

Brief description of the drawings.

In fig. Figure 1 shows a general view of the system for automatic search and correction of the object’s rotation axis.

In fig. 1 is indicated: 1 - X-ray source, 2 - object under study, 3 - device for rotating and positioning the object under study (the dotted line indicates the axis of rotation), 4 - flat panel detector, 5 - computing device for collecting and processing X-ray images.

In fig. Figure 2 shows a block diagram of data processing.

In fig. Figure 3 shows a block diagram of a two-stage determination of the angle of inclination of the rotation axis.

In fig. Figure 4 shows a block diagram of a two-stage determination of the amount of axis shift during rotation of the average projection image.

In fig. Figure 5 shows a graph of the function for assessing the asymmetry of the rotated preprocessed averaged projection image l(Δx) for various trial values of the shift in the position of the Δx axis for the averaged projection image of the phantom.

In fig. Figure 6 shows a graph of the function for assessing the asymmetry of the rotated preprocessed averaged projection image l(Δx) for various trial values of the shift in the position of the Δx axis for the averaged projection image of the phantom in the vicinity of the minimum point.

In fig. Figure 7 shows a graph of the function for assessing the asymmetry of the rotated preprocessed averaged projection image l(Δx) for various trial values of the shift in the position of the Δx axis for the averaged projection image of a tube with chess.

In fig. Figure 8 shows a graph of the function for assessing the asymmetry of the rotated preprocessed averaged projection image l(Δx) for various trial values of the shift in the position of the Δx axis for the averaged projection image of a tube with chess in the vicinity of the minimum point.

In fig. Figure 9 shows a block diagram for refining the magnitude of the shift of the rotation axis based on the initial averaged projection image.

In fig. Figure 10 shows a graph of the objective functional for synthetic data, calculated from an unscaled average projection image.

In fig. Figure 11 shows a graph of the objective functional for real data, calculated from an unscaled average projection image.

In fig. Figure 12 shows the reconstructed central layer of synthetic data (phantom with geometric figures): a.) ideal image, b.) reconstruction of the central layer without correction of the rotation axis, c.) reconstruction of the central layer with the corrected rotation axis.

In fig. Figure 13 shows the reconstructed central layer of real data (tube with chessboard): a.) ideal, reconstruction of the central layer with manually selected parameters for the position of the rotation axis, b.) reconstruction of the central layer without correction of the rotation axis, c.) reconstruction of the central layer with the corrected rotation axis .

The operating diagram of the claimed system is presented in Fig. 1.

The system consists of an X-ray source 1, a flat-panel detector 4, between which a device for rotating and positioning 3 of the object under study 2 and a computing device 5 connected to the detector 4 are installed.

The scanning pattern is circular. The detector is flat panel. The geometric model of radiation propagation is parallel, fan-shaped or conical, the detector dimensions are from 100×100 to 5000×5000 pixels, the angular size of the object is up to 30 degrees. Source - synchrotron or laboratory.

During scanning, the data obtained by the flat panel detector is transferred to a computing device, which carries out the following data processing steps:

- obtaining projection data,

- calculation of the average projection image,

- preprocessing of the average projection image using linear normalization and scaling,

- two-stage search for the values of the angle of inclination and shift of the axis of rotation,

- clarification of the position of the object’s rotation axis.

Receiving Projection Data

In the first stage, the computing device receives projection data from the flat panel detector.

Calculation of the average projection image.

Calculation of the arithmetically averaged projection image :

where p _ϕ =p _ϕ (x, z) is a projection image of an object (single-channel image) corresponding to the projection angle ϕ ∈ Ф, Ф is the set of viewed projection angles in the amount of |Ф|.

Let us assume that the averaged projection image is an axisymmetric image, and the axis of symmetry coincides with the projection of the object's rotation axis onto the plane of the detector window when illuminated by an X-ray source. In the case of shifts, the averaged projection image can be considered an approximately axisymmetric image: the asymmetry of the averaged projection image becomes noticeable at shifts exceeding several tens to hundreds of detector pixels. In a parallel geometric scheme, when considering the set of permissible movements, the axis of rotation of the averaged projection image is always strictly axisymmetric.

The average projection image consists of a set of ellipses - projections of circular trajectories of movement of individual material points of the object under study during rotation around an axis. The latter serves as a justification for the presence of axial symmetry of the average projection image.

Under the assumption of axisymmetricity, the symmetry axis of the averaged projection image coincides with the projection of the object’s rotation axis onto the plane of the detector window: the angle of inclination of the axis in a plane parallel to the plane of the detector window is equal to the angle between the symmetry axis of the averaged projection image and the vertical, and the shift of the rotation axis differs from the shift the symmetry axis of the average projection image from the center of the image by the number of times specified by the magnification factor of the conical geometric scheme, or coincides with the shift of the symmetry axis of the average projection image in the parallel geometric scheme. Thus, the problem of determining the scalar parameters of the position of the object’s rotation axis - the magnitude of the shift and tilt in a plane parallel to the plane of the detector window - is reduced to the problem of finding the equation of the symmetry axis of the average projection image, in particular, its tilt angle and shift in the horizontal direction (from the center Images).

Preprocessing the average projection image using linear normalization and scaling.

After calculating the average projection image it is preprocessed. The image preprocessing stage includes linear normalization and scaling.

First, linear minimax normalization of the average projection image is carried out, in which the pixel intensity values of the average projection image are reduced to the [0,1] scale.

The average projection image is then scaled in the horizontal and vertical directions by the same scaling factor When scaling, the dimensions of the average projection image are reduced so that their ratio remains constant, and the smallest dimension of the scaled (with averaging) image is 128 pixels. As part of the scaling step of the average projection image, interpolation is performed.

Scaling is used to speed up subsequent data processing steps, which include enumerating different rotation axis shift values across the entire width of the image.

We will distinguish between preprocessed and unpreprocessed averaged projection images and use different notations for these objects: original unpreprocessed image and preprocessed (normalized and scaled) image

Two-stage search for the values of the angle of inclination and shift of the axis of rotation

After preprocessing the averaged projection image, a two-stage search is performed using the enumeration method over a uniform sparse grid and a uniform grid with a reduced step of the values of the inclination angle of the rotation axis in a plane parallel to the plane of the detector window (all considered grids of parameter values, unless otherwise stated, are considered uniform; instead of uniform grids, computational grids of some other classes can be used). While estimating the axis tilt angle, an approximation to the desired value of the rotation axis shift is also found, which is refined in the next step.

Stage 1 consists of searching for the angle of inclination of the rotation axis using a sparse grid of values; within stage 2, the value of the inclination angle of the axis found at stage 1 is refined using a grid with a reduced step.

Let us describe in detail each of the two stages of the current step of the method. The steps are visualized through the flowchart shown in FIG. 3.

Stage 1: searching for the inclination angle of the rotation axis using a sparse grid of values.

Rotate the preprocessed average projection image. We have a preprocessed average projection image We will look for the first approximation to the desired value of the angle of inclination of the rotation axis in the interval from start_angle to end_angle degrees with a step of first_step.

The parameters first_step, start_angle and end_angle of the sparse uniform grid of values under consideration can be chosen to be 1.0, -15.0 and 15.0. To search for a first approximation of the axis tilt angle, the grid step is chosen not to exceed 1 degree.

For each angle value from the considered interval, the preprocessed image is rotated by angle degrees. Rotation is carried out relative to the center of the image using interpolation.

The preprocessed image and the original averaged projection image rotated by angle degrees will be denoted by And respectively. In such notation, in particular,

Note that due to symmetry relative to the projection of the rotation axis onto the plane of the detector window, if angle is the correctly found value of the axis inclination angle, has an exactly vertical axis of symmetry. Consequently, the best approximation angle of the true value of the axis inclination angle corresponds to the situation when the function for estimating the vertical asymmetry of the image, defined below minimal.

Target functionality. Function for estimating the vertical asymmetry of the average projection image. To measure the function for assessing the vertical asymmetry of an arbitrary image P, we define l(Δx) - the function for assessing the asymmetry of the average projection image P relative to the vertical, shifted by Δx pixels relative to the middle vertical of the image in the horizontal direction (target functional):

where F is the operator of mirror image reflection in the horizontal direction, S _Δx is the operator of shifting the image by Δx in the horizontal direction, ° is the operation of composition of images, l _q is the standard q Minkowski norm. Instead of the Minkowski norm, another matrix norm can be calculated.

The vertical axis, shifted horizontally relative to the center of the image by Δx pixels, is the axis of symmetry of the image P if, with horizontal mirroring of P and a subsequent shift by 2⋅Δx pixels, we also obtain the original image P. The evaluation function l(Δx) of the image asymmetry relative to the axis shifted from center by Δx pixels, there is l _q - the norm of the difference between the original image P and its mirrored copy shifted by 2⋅Δx pixels (S _2⋅Δx °F)[P]. The smaller the value of l(Δx), the more confident we can be that x=Δx is the vertical axis of symmetry of the image (counting from the center of the image).

The function for assessing vertical asymmetry l _P of image P shows how likely image P is to have any vertical axis of symmetry:

where the minimum is taken over some computational grid possible values of shifts from the center P of the assumed vertical axis of symmetry.

The value l _P =0 (note, l _P ≥0) of the function for assessing vertical asymmetry indicates the existence of a vertical axis of symmetry in the average projection image P. The greater the value of 1 _P , the less symmetrical the image (assuming - a fairly fine mesh).

We consider the vertical axis of symmetry of the average projection image P to be such a vertical axis shifted relative to the center of the image by Δx pixels, on which the functional l _P (Δx) reaches the value l _P :l _P (Δx)=l _P - the condition that the vertical axis in question is the vertical axis of symmetry of the averaged projection image R.

Two-stage minimization of the objective functional. The desired value of the angle of inclination of the rotation axis is such that the rotated preprocessed average projection image has a vertical axis of symmetry, that is, the optimal angle value is selected from the condition of minimizing the function for assessing vertical asymmetry average projection image

For each angle value of the sparse mesh (from start_angle to end_angle degrees with step first_step) the value is calculated during a two-stage minimization of the objective functional (Δx) (see Fig. 4).

Stage 1: minimizing the target functional over a sparse grid of values Δx from -width to width pixels with a large step first_step, where width is the width of the image , and the first step can be set to 5 (at least a few pixels). At the end of stage 1, we obtain an initial approximation of the minimum value of the vertical asymmetry estimation function and the value Δх, which specifies the position of the minimizing vertical axis of symmetry. Graphs of target functionality (Δх) (at the angle value coinciding with the true value of the angle of inclination of the rotation axis) are shown in Fig. 5, 6, 7, 8.

Stage 2: minimization of the target functional using a refined grid of shift values Δx, during which the target functional is minimized sequentially for each p ∈ {1,0, -1} (Δх) in the vicinity of the previously found minimum point Δх=best_shift of radius 10 ^р with a step of 10 ^р-1 (the step is equal to the linear pixel size at р=1, a tenth and a hundredth of the linear pixel size at р=0 and р=-1, respectively, minimization by direct search). Examples of graphs of the objective functional in the vicinity of the minimum point can be seen in Fig. 6 and 8.

The result of stage 2 is the calculated value and Δх, at which the target functional (Δх) is equal (when the preprocessed image is previously rotated by an angle relative to the center of the average projection image).

Next, we rotate the preprocessed averaged projection image to another angle (see Fig. 3), and calculate the function for estimating vertical asymmetry for the generated rotation angle angle If the found norm value of the vertical asymmetry estimation function of the rotated averaged projection image is less than the estimation function for the previous rotation (that is, min_norm > norm), we update the values of the tilt angle and axis shift angle_best=angle, shift_best=shift and min_norm=norm. The values of the angle of rotation are sorted from start_angle to end_angle degrees with a step of first_step.

At this stage, stage 1 of the two-stage search for the values of the angle of inclination and shift of the rotation axis is completed. At the end of stage 1, we have initial approximations of the angle of inclination of the rotation axis and the shift of the rotation axis, at which the corresponding optimal value of the vertical asymmetry estimation function is achieved.

Stage 2: clarifying the tilt value of the rotation axis.

As part of stage 2 of the step of calculating the value of the angle of inclination of the rotation axis, refined grids of values of the inclination angle in the vicinity of the optimal value angle_best found during stage 1 are considered. Using refined meshes, the function of estimating the vertical asymmetry of the rotated preprocessed averaged projection image is minimized

Let's describe stage 2 in more detail. For each p ∈ {0, -1}, a grid of angle values is considered from angle best -2⋅10 ^p to angle_best +2⋅10 ^p with a step of 10 ^p-1 (a tenth of a degree for p=0 and a hundredth of a degree for p= -1), where angle_best is the previously determined optimal value of the axis tilt angle. For each angle, the vertical asymmetry evaluation function is calculated through the already described two-stage minimization of the target functional. If the new value of the vertical asymmetry estimation function is less than previously obtained for other angle values, we update the values of angle_best=angle, shift_best=shift and min_norm=norm.

The result of stage 2 is the optimal value of the angle of inclination of the rotation axis angle_best, refined using refined meshes. For the calculated angle_best value, the value of the shift of the rotation axis shift_best is also found, at which the found minimum value of the function for assessing the asymmetry of the average projection image is achieved. The value angle_best is the answer to the question about the value of the angle of inclination of the rotation axis and is not further specified.

The shift_best value serves as an initial approximation of the desired value of the rotation axis shift, which is then subject to refinement.

Clarifying the magnitude of the shift of the rotation axis using an averaged projection image.

In the previous step, the values of the inclination angle of the rotation axis and the shift of the rotation axis were calculated. At the step of clarifying the value of the shift of the rotation axis, the previously found shift_best value is refined.

The axis shift value is refined using the original unpreprocessed average projection image (see Fig. 9).

As part of the step of clarifying the value of the shift of the rotation axis using the original averaged projection image, the steps outlined below are performed.

1) Linear minimax normalization of the image is performed (pixel intensity values of the average projection image are scaled to [0,1]). For uniformity, we denote the result of normalization of the original averaged projection image We emphasize that the average projection image is not scaled.

2) The normalized average projection image is rotated around the center of the image by the found angle of inclination of the rotation axis, equal to angle_best. We will denote the result of the rotation The rotated normalized averaged projection image has a vertical axis of symmetry, which we will look for in the vicinity of the position of the vertical axis Since the initial approximation of the axis shift value was obtained using a scaled average projection image, shift_best is divided by the scaling_factor.

3) Consider the grid of possible values of the shift of the axis of rotation from before pixels with a small pitch equal to a hundredth of the linear pixel size. Using this grid, the value of the target functional is minimized by direct enumeration of the shift values Δx The minimum point shift_best is the final answer of the method to the question about the value of the shift of the rotation axis. Examples of graphs of the target functional, the values of which are calculated for various values of the shift of the rotation axis Δx using a normalized unscaled image, are shown in Fig. 10 and 11.

After calculating shift_best, the specifics of the X-ray propagation model should be taken into account, namely, in the cone model, it is necessary to divide shift_best by the cone magnification factor (the same applies to the fan-shaped scheme). The parallel model of X-ray propagation does not require normalization of the shift value of the object's rotation axis found from the average projection image.

So, the values of the angle of inclination and shift of the rotation axis have been found.

Rotation axis position correction

After finding the parameters of the position of the rotation axis (the values of the angle of inclination and shift of the rotation axis), each projection image r _f , ϕ ∈ Ф, is rotated relative to the center of the image by the angle of inclination of the rotation axis angle_best. Then each rotated projection image is shifted by shift_best pixels. When calculating rotations and shifts of projection images, interpolation is used.

Next, a tomographic reconstruction algorithm can be applied to obtain an image of the internal structure of the object under study.

An example of the presented method

The results of the invention on synthetic and real data are presented in Fig. 12 and 13. The figures show reconstructions of the central tomographic layers of synthetic data (phantom with geometric figures) and real data (tube with chessboard) with and without correction of the rotation axis by the proposed automatic method for aligning the rotation axis (Figs. 12b, 12c and 13b, 13c ). In fig. 12a, 13a show ideal images of the central layer of the phantom and the tube with chessboard, with which the resulting reconstructed central layers are compared. To measure changes in the quality of the resulting reconstructions, the L ₂ norm of the difference between the reconstruction and the ideal image was calculated: if the reconstructed image is represented by a matrix and the ideal image is a matrix then the value l ₂ of the image difference norm is equal to The Feldkamp-Davis-Kress algorithm was used for reconstruction. The value of l ₂ norms for the difference between reconstructions and ideal images is presented in Table 1.

Synthetic projection images (packet of 1024 × 1024 × 720 voxels) are simulated with a rotation axis tilt angle of 9.20 degrees in a plane parallel to the plane of the detector window, and a shift of 2.83 pixels. The axis tilt angle value found by the proposed method is 9.20 degrees, the axis shift value found is 2.81 pixels. The result of correcting the position of the rotation axis is shown in Fig. 12th century The reconstruction of the central layer when correcting the rotation axis is visually similar to the ideal image of the central layer of the phantom (cf. 12a and 12c). The lower value of L ₂ norm of the difference between the reconstruction and the ideal image in the case of using rotation axis correction confirms the improvement in the quality of reconstruction when using the proposed axis correction tool (Table 1, cf. Figs. 12b and 12c).

When applied to a package of real data of a tube with a chessboard (package of 2064 × 1548 × 1024 voxels), an improvement in the quality of reconstruction is visually obvious (Fig. 13b and 13c). The found value of the inclination angle of the rotation axis is 0.33 degrees, the found value of the axis shift is -7.09. Manual selection of the rotation axis position parameters confirms the accuracy of the method results: the optimal shift value is between -7.00 and -8.00 pixels, while the optimal tilt angle is between 0.00 and 0.50 degrees. In fig. Figure 13a shows the reconstruction of the central layer with optimal axis parameters, selected manually and equal to -7.00 pixels and 0.30 degrees.

The examples given show that the use of the described means of automatically determining the position of the rotation axis makes it possible to obtain tomographic reconstructions that are closer to the ideal image, both in visual perception and according to the measured L ₂ norm (see Table 1), than in the case without axis correction.

The use of the proposed invention precedes the step of reconstructing the internal structure of the object. Correcting the position of the rotation axis before subsequent reconstruction makes it possible to obtain a reconstructed image of higher quality.

In this case, the described invention can also be used as a means of preliminary calibration of measuring equipment using special phantoms. This application of the invention makes it possible to improve the quality of results in subsequent target measurements.

Claims (6)

A system for automatically searching and correcting the axis of rotation of an object when performing X-ray tomography, containing an X-ray source, a device for rotating and positioning the object, a flat panel detector and a computing device, characterized in that the flat panel detector is configured to generate projection data of the object, and the computing device is configured to : a) obtaining projection data, b) calculating the average projection image, c) preprocessing the average projection image using linear normalization and scaling, d) two-stage search for the values of the angle of inclination and shift of the axis of rotation, e) clarifying the position of the object’s rotation axis.