CN103211608B - Computational methods based on the spiral CT image reconstruction mesoporous weighting that GPU accelerates - Google Patents

Abstract

The present invention provides the computational methods of a kind of spiral CT image reconstruction mesoporous weighting accelerated based on GPU, comprises the steps: that a. determines the aperture weighted value of all asymmetric radiation of each voxel；B. the aperture weighted value of asymmetric radiation all to each voxelIt is added, obtains the aperture weighted value of this voxel；C. the aperture weighted value of all voxels is normalized calculating, obtains the aperture weighted value after each voxel normalization.The present invention effectively solves the uneven slip of edge detector unit detection data by the method that aperture weights, and need not conditional statement and just can accurately calculate multislice spiral CT image reconstruct mesoporous weighting and normalized value, thus improve quality and the efficiency of restructing algorithm of reconstruct image.

Description

Computational methods based on the spiral CT image reconstruction mesoporous weighting that GPU accelerates

Technical field

The present invention relates to Medical Image Processing, particularly relate to the computational methods of a kind of spiral CT image reconstruction mesoporous weighting accelerated based on GPU.

Background technology

The image reconstruction speed of multi-layer spiral CT is always a focus of CT technical research, and technology relatively common in existing accelerated method can be divided into three types:

The first type is accelerated calculating, as disclosed a kind of fast parallel cone-beam reconstruction system and method using one or more microprocessor CPU in the Chinese patent of Patent No. CN02805089.4 for utilizing CPU cluster.But the method is just for cone beam reconstruction, does not mention the reconstruction model that multi-layer spiral CT is relevant, and CPU is expensive, will reach reasonable acceleration effect cost the most too high.

The second type is accelerated for utilizing FPGA this kind of specific calculations hardware device, as: the method disclosed in article " HighSpeedCTImageReconstructionusingFPGA " (InternationalJournalofComputerApplications (0,975 8887) Volume22 No.4, May2011).But the method the most simply mentions basic CT reconstruct, not relating to multi-layer spiral CT correlation technique, and FPGA this specific calculations hardware transplantability after algorithm designs is not strong, the dependence to hardware is serious.

The third type is accelerated for utilizing graphic process unit (GPU).The cost of this mode reduces a lot for CPU cluster, and the floating-point disposal ability of existing GPU is far above CPU, and only one piece of video card that it needs, transplantability and compatible specific calculations hardware device more this than FPGA are the best, and the acceleration carrying out CT restructing algorithm hence with GPU can reach faster reconstructed velocity under conditions of low cost.Also compare many for the patent in terms of this, as the Chinese patent of Patent No. CN200810114478.1 discloses a kind of CT parallel reconstructing system and formation method, utilize GPU cluster to carry out CT reconstruct；The Chinese patent of Patent No. CN200810113846.0 discloses the GPU accelerated method of a kind of CT image reconstruction, but both approaches is not particularly suited for multislice spiral CT image reconstruct, does not the most comprise the computational problem to aperture weighting function.Additionally, although the Chinese patent of Patent No. CN200910248774.5 discloses one and utilizes graphic process unit, it is achieved the method for 3 D back projection, the method is applicable to multislice spiral CT image reconstruct, but the calculating of the most not mentioned aperture weighted value, and normalized computational problem.

In multislice spiral CT image reconstructs, the quality calculating the efficiency on restructing algorithm and reconstruct image of aperture weight values suffers from very important impact.During Spiral CT scan, patient moves all along being parallel to detector direction, can there is the mutation process that detectable signal grows out of nothing for several detector units at detector edge, and then impact detects the concordance of data.It is thus desirable to add an aperture weighting function to smooth this process grown out of nothing, thus reducing arch artifact, the importance calculated aperture weighting function in article " WeightedFBP asimpleapproximate3DFBPalgorithmformulti-slicespiralCTwi thgooddoseusageforarbitrarypitch " (Phys.Med.Biol.49 (2004) 2,209 2218) has detailed description.But this article is the theory analysis that this weighting function carries out basis, does not relate to how to implement, does not more mention and how to utilize GPU to be accelerated.

Therefore, the necessary computational methods that a kind of multislice spiral CT image reconstruct mesoporous weighting accelerated in graphic process unit (GPU) is provided, can solve the problem that the uneven slip of edge detector unit detection data, and need not conditional statement and just can accurately calculate multislice spiral CT image reconstruct mesoporous weighting and normalized value, improve quality and the efficiency of restructing algorithm of reconstruct image.

Summary of the invention

It is an object of the invention to provide the computational methods of a kind of spiral CT image reconstruction mesoporous weighting accelerated based on GPU, can effectively solve the problem that the uneven slip of edge detector unit detection data, and need not conditional statement and just can accurately calculate multislice spiral CT image reconstruct mesoporous weighting and normalized value, thus improve quality and the efficiency of restructing algorithm of reconstruct image.

For achieving the above object, the present invention adopts the following technical scheme that the computational methods of a kind of spiral CT image reconstruction mesoporous weighting accelerated based on GPU, comprises the steps: that a. determines the aperture weighted value of all asymmetric radiation of each voxel；B. the aperture weighted value of asymmetric radiation all to each voxelIt is added, obtains the aperture weighted value of this voxel；C. the aperture weighted value of all voxels is normalized calculating, obtains the aperture weighted value after each voxel normalization。

Further, described aperture weighted valueIt is calculated as follows；

(1)

Wherein, Q is the parameter of controlling curve smoothness, and q is the detector coordinate in Z-direction,For current projection angle, N is the row of detector,For the voxel after reconstruct at X, Y, the coordinate of Z-direction.

Further, described detector is calculated by formula (2) at the coordinate q of Z-direction:

(2)

Wherein,For current projection angleUnder scanning bed physical location in a z-direction, t is the coordinate in probe access direction:, R is the distance of x-ray source and center of rotation,For the scanning bed distance of movement in a circle scanning,,For center of rotation in the value of port number corresponding to channel direction projection.

Further, described spiral CT is multi-layer spiral CT, and the asymmetric radiation differing 180 ° or 360 ° in helical scanning, when the coordinate q of Z-direction, is divided into the asymmetric radiation of odd-numbered and the asymmetric radiation of even-numbered, and is respectively calculated by calculating detector.

Further, described even-numbered is the asymmetric radiation of the n coordinate in detector z direction(3) calculate as follows, and wherein n is the integer more than or equal to 0；

(3)

Initial value(4) calculate as follows:

(4)

Being the spacing between two neighbouring even-numbered numbering asymmetric radiation, (5) calculate as follows:

(5).

Further, described odd-numbered is the asymmetric radiation of the n coordinate in detector z direction(6) calculate as follows, and wherein n is the integer more than or equal to 0；

(6)

Initial value(7) are calculated as follows:

(7)

Being the spacing between two adjacent odd numbering asymmetric radiation, (8) are calculated as follows:

(8).

Further, described even-numbered n scope is。

Further, described odd-numbered n scope is。

Further, the aperture weighted value after described normalization(11) calculate as follows:

(11)

Wherein, i is positive integer, and 2n is the quantity through all asymmetric radiation of this voxel.

Further, the aperture weighted value after described normalizationIt is stored in the shared video memory of GPU.

The present invention contrasts prior art following beneficial effect: the computational methods of the spiral CT image reconstruction mesoporous weighting accelerated based on GPU that the present invention provides, the method weighted by aperture effectively solves the uneven slip of edge detector unit detection data.In addition, the asymmetric radiation differing 180 ° or 360 ° in multi-layer spiral CT helical scanning is divided into the asymmetric radiation of odd-numbered and the asymmetric radiation separate computations of even-numbered by the present invention, improve the efficiency of restructing algorithm further, solve the high consumption characteristics that conditional statement is processed by GPU.

Accompanying drawing explanation

Fig. 1 is the flow chart of steps of the computational methods of the spiral CT image reconstruction mesoporous weighting that the present invention accelerates based on GPU.

Fig. 2 is the relation schematic diagram of the aperture weighted value in spiral CT image reconstruction of the present invention and the coordinate of detector Z-direction.

Fig. 3 is the arrangement schematic diagram of asymmetric radiation in multi-slice Spiral CT of the present invention.

Fig. 4 is the multi-layer spiral CT reformatted slices figure to actual die body in the embodiment of the present invention.

Detailed description of the invention

The invention will be further described with embodiment below in conjunction with the accompanying drawings.

Fig. 1 is the flow chart of steps of the computational methods of multislice spiral CT image of the present invention reconstruct mesoporous weighting, and Fig. 2 is the relation schematic diagram of the aperture weighted value in multislice spiral CT image of the present invention reconstruct and the coordinate of detector Z-direction.

Refer to Fig. 1, the computational methods of the multislice spiral CT image reconstruct mesoporous weighting accelerated based on GPU that the present invention provides, comprise the steps:

First, in step S101, (1) determines the aperture weighted value of all asymmetric radiation of each tissue points as follows；

(1)

Aperture weighted valueFor being multiplied by a coefficient on each probe unit of detector Z-direction, this coefficient is less than or equal to 1, and the value of central passage is 1, and the value of edge gateway is the numerical value less than 1.In above-mentioned formula (1), Q is the parameter of controlling curve smoothness, and q is the detector coordinate in Z-direction,For current projection angle, N is the number of plies of the row of detector, i.e. multi-layer spiral CT,For the i-th voxel after reconstruct at X, Y, the coordinate of Z-direction；Detector is in coordinate q and the aperture weighted value of Z-directionRelation as shown in Figure 2.

In Three-dimensional Multi-slice Spiral CT image reconstruction, after data rearrangement, described detector can (2) be calculated as follows at the coordinate q of Z-direction:

(2)

In above-mentioned formula (2),For i-th voxel at the coordinate of z-axis,For current projection angleUnder scanning bed physical location in a z-direction, t is the coordinate in probe access direction:, R is the distance of x-ray source and center of rotation,For the scanning bed distance of movement in a circle scanning,,For center of rotation in the value of port number corresponding to channel direction projection.

Owing to there is a plurality of asymmetric radiation in helical scanning, and the asymmetric radiation of difference 180 ° and 360 ° all there are differences in spacing and initial position, so asymmetric radiation is divided into the asymmetric radiation of odd-numbered and the asymmetric radiation of even-numbered, and is respectively calculated.

Described even-numbered is the asymmetric radiation of the n coordinate in detector z direction(3) calculate as follows, and wherein n is the integer more than or equal to 0；

(3)

Initial value(4) calculate as follows:

(4)

Being the spacing between two neighbouring even-numbered numbering asymmetric radiation, (5) calculate as follows:

(5).

Similarly, described odd-numbered is the asymmetric radiation of the n coordinate in detector z direction(6) calculate as follows, and wherein n is the integer more than or equal to 0:

(6)

Initial value(7) are calculated as follows:

(7)

Being the spacing between two adjacent odd numbering asymmetric radiation, (8) are calculated as follows:

(8).

In order to determine the scope of odd, even number numbering n, also carry out according to the mode of odd even calculated hole diameters weights when, first analyze the situation of even-numbered, it is assumed that the initial value position of even-numbered is No. 0 asymmetric radiation, as shown in Figure 3.

Owing to detector origin coordinates is 0, termination coordinate be N, N be the number of plies of multi-layer spiral CT, therefore can be calculated as follows according to these parameters and obtain detector original positionPositional value in-2 to+2 these sequences:

(9)

With detector final positionPositional value in this sequence:

(10).

Thus can accurately calculate the numbering of each asymmetric radiation,For being not less thanSmallest positive integral,For being not more thanMaximum integer, just can accurately calculate each even number asymmetric radiation in formula (3) after obtaining this valueValue, it is not necessary to any judgement operation, it is only necessary to twice plus/minus and twice division.

In like manner, odd-numbered asymmetric radiation can be calculated, after obtaining this value, just can accurately calculate each odd-numbered symmetry ray in formula (6)Value.

Then, in step s 102, the aperture weighted value to all asymmetric radiation of each voxel calculated in step S101It is added, obtains the aperture weighted value of this voxel；

Finally, in step s 103, the aperture weighted value to all tissue points obtained in step S102, (11) are normalized calculating as follows, obtain the aperture weighted value after each tissue points normalization:

(11)

Wherein, i is positive integer, and 2n is the quantity through all asymmetric radiation of this voxel.

Due to these aperture weightsValue for the probe unit of all Z-directions is all identical, therefore the number of plies according to detector Z-direction when program starts is had only to, calculate once, conduct interviews owing to each GPU thread is required for transferring this one piece of data, it is contemplated that to the different types of video memory of GPU, the difference of its access efficiency, the present embodiment is by calculated aperture weightsIt is stored in the shared video memory of GPU.Shared video memory is positioned in the core processor of each GPU, belongs to on-chip memory cell, and its access speed is the fastest in a few class GPU video memory.

To illustrate and verify the effect of the present invention with concrete data and reconstructed image quality below.

The data obtained with the scanning of actual die body carry out 3-D view reconstruct and obtain 64 slice maps, take a slice map therein to contrast, one of which data are not for having calculated hole diameters weighting only to carry out normalized, and another set data calculate aperture and weight and be normalized.The parameters that this slice map uses is:, 64 Layer Detection devices, detector is in spacing 0.625mm in z direction, and reconstruct visual field size is 210mm, and light source to center of rotation distance is 570mm, and reconstruct thickness is 0.625, and Q-value is 0.8.

Fig. 4 is the multi-layer spiral CT reformatted slices figure to actual die body, and Fig. 4 a is the reformatted slices figure obtained under not having calculated hole diameters weighting conditions, and Fig. 4 b is the reformatted slices figure obtained under calculated hole diameters weighting conditions.Fig. 4 c is the differential chart of Fig. 4 a and Fig. 4 b.

From fig. 4, it can be seen that there is a little arch artifact in the upper left corner of Fig. 4 a, can be evident that in Fig. 4 c the slice map that two distinct methods obtain has significant diversity.Artifact inconspicuous in Fig. 4 a, this is owing in the data actually obtained, noise level is higher, masking partially due to do not add the arch artifact that aperture weighted calculation is brought, the artifact that this part is blanked is in the differential chart 4c of two methods or can clearly reflect.This most just illustrates that the computational methods that the spiral CT image reconstruction mesoporous based on GPU acceleration that the present invention provides weights have critically important impact for the picture quality that actual multi-layer spiral CT reconstructs.

Underneath with concrete data, the computational efficiency that the multi-layer spiral CT reconstruct mesoporous based on GPU acceleration that the present invention provides weights is described.

Experiment condition is identical with the parameter used in above-mentioned Fig. 4 reformatted slices figure, the experiment condition of contrast groups and the GPU code that use are completely the same with the computational methods that the present invention provides, and only difference is that in contrast groups that calculated hole diameters weights have employed and judge that statement calculates.Shown in the following form of computational efficiency:

Reconstruct number of voxel	512x512x1	512x512x16	512x512x32	512x512x64
					Contrast groups (millisecond)	796	5425	10847	21719
The inventive method (millisecond)	677	3081	6153	12378

Can clearly find out that the computational methods that the spiral CT image reconstruction mesoporous based on GPU acceleration that the present invention provides weights improve by about one time in computational efficiency from upper table, this is to have given up completely due to the computational methods of present invention offer to judge that statement can accurately calculate quantity and the value of its aperture weighting of the asymmetric radiation required for the tissue points of each reconstruct.

To sum up, the computational methods of the spiral CT image reconstruction mesoporous weighting accelerated based on GPU that the present invention provides, the method weighted by aperture effectively solves the uneven slip of edge detector unit detection data.In addition, the asymmetric radiation differing 180 ° or 360 ° in multi-layer spiral CT helical scanning is divided into the asymmetric radiation of odd-numbered and the asymmetric radiation separate computations of even-numbered by the present invention, improve the efficiency of restructing algorithm further, solve the high consumption characteristics that conditional statement is processed by GPU.

Although the present invention discloses as above with preferred embodiment; so it is not limited to the present invention, any those skilled in the art, without departing from the spirit and scope of the present invention; when making a little amendment and perfect, therefore protection scope of the present invention is when with being as the criterion that claims are defined.

Claims (7)

1. the computational methods of the spiral CT image reconstruction mesoporous weighting accelerated based on GPU, it is characterised in that comprise the steps:

A. the aperture weighted value W (θ of the asymmetric radiation of each voxel is determined_j, q (x_i, y_i, z_i)), described θ_jFor current projection angle, described x_i, y_i, z_iFor the i-th voxel after reconstruct at X, Y, the coordinate of Z-direction；

Described aperture weighted value W (θ_j, q (x_i, y_i, z_i)) it is calculated as follows:

Wherein, Q is the parameter of controlling curve smoothness；

B. the aperture weighted value W (θ of asymmetric radiation all to each voxel_j, q (x_i, y_i, z_i)) be added, obtain the aperture weighted value of this voxel；

C. the aperture weighted value of all voxels is normalized calculating, obtains the aperture weighted value after each voxel normalizationDescribed q is the detector coordinate in Z-direction, and N is the row of detector；

Described detector is calculated by formula (2) at the coordinate q of Z-direction:

$q=\frac{\left(z_i-{\mathrm{bedpos}}_j-\arcsin \left(t/R\right)\frac{pitch\_bed}{2\mathrm{\π}}\right)R}{\sqrt{R^2-t{(x_i,y_i,{\mathrm{\θ}}_j)}^2}+v(x_i,y_i,{\mathrm{\θ}}_j)}+{\mathrm{Center}}_{Channel}---\left(2\right)$

Wherein, bedpos_jFor current projection angle θ_jUnder scanning bed physical location in z-direction, t is the coordinate in probe access direction: t=y_icosθ_j-x_isinθ_j, R is the distance of x-ray source and center of rotation, and pitch_bed is the scanning bed distance of movement, v=x in a circle scanning_icosθ_j+y_isinθ_j, Center_channelFor center of rotation in the value of port number corresponding to channel direction projection.

2. the computational methods of the spiral CT image reconstruction mesoporous weighting accelerated based on GPU as claimed in claim 1, it is characterized in that, described spiral CT is multi-layer spiral CT, calculating detector is when the coordinate q of Z-direction, the asymmetric radiation differing 180 ° and 360 ° in helical scanning is divided into the asymmetric radiation of odd-numbered and the asymmetric radiation of even-numbered, and is respectively calculated.

3. the computational methods of the spiral CT image reconstruction mesoporous weighting accelerated based on GPU as claimed in claim 2, it is characterised in that described even-numbered is the asymmetric radiation coordinate in detector Z-direction of n(3) calculate as follows, and wherein n is the integer more than or equal to 0；

$q_n^{even}=q_{start}^{even}+n\×q_{\mathrm{int}erval}^{even}---\left(3\right)$

Initial value(4) calculate as follows:

$q_{start}^{even}=\frac{(z_i-{\mathrm{bedpos}}_j-\arcsin (t/R\left)\right)R}{\sqrt{R^2-t{(x_i,y_i,{\mathrm{\θ}}_j)}^2}+v(x_i,y_i,{\mathrm{\θ}}_j)}+{\mathrm{Center}}_{Channel}---\left(4\right)$

Being the spacing between two neighbouring even-numbered numbering asymmetric radiation, (5) calculate as follows:

$q_{\mathrm{int}erval}^{even}=\frac{-pitch\_bed\×R}{\sqrt{R^2-t{\left(x_i,y_i,{\mathrm{\θ}}_j\right)}^2}+v\left(x_i,y_i,{\mathrm{\θ}}_j\right)}---\left(5\right).$

4. the computational methods of the spiral CT image reconstruction mesoporous weighting accelerated based on GPU as claimed in claim 3, it is characterized in that, described odd-numbered is that the asymmetric radiation of n calculates at the coordinate (6) as follows of detector Z-direction, and wherein n is the integer more than or equal to 0；

$q_n^{odd}=q_{start}^{odd}+n\×q_{\mathrm{int}erval}^{odd}---\left(6\right)$

Initial value(7) are calculated as follows:

$q_{start}^{odd}=\frac{(z_i-{\mathrm{bedpos}}_j-pitch\_bed/2-\arcsin (t/R\left)\right)R}{\sqrt{R^2-t{(x_i,y_i,{\mathrm{\θ}}_j)}^2}-v(x_i,y_i,{\mathrm{\θ}}_j)}+{\mathrm{Center}}_{Channel}---\left(7\right)$

Being the spacing between two adjacent odd numbering asymmetric radiation, (8) are calculated as follows:

$q_{\mathrm{int}erval}^{odd}=\frac{-pitch\_bed\×R}{\sqrt{R^2-t{(x_i,y_i,{\mathrm{\θ}}_j)}^2}-v(x_i,y_i,{\mathrm{\θ}}_j)}---\left(8\right).$

5. the computational methods of the spiral CT image reconstruction mesoporous weighting accelerated based on GPU as claimed in claim 4, it is characterised in that described even-numbered n scope isWhereinFor being not less thanSmallest positive integral,For being not more thanMaximum integer.

6. the computational methods of the spiral CT image reconstruction mesoporous weighting accelerated based on GPU as claimed in claim 4, it is characterised in that described odd-numbered n scope isWhereinFor being not less thanSmallest positive integral,For being not more thanMaximum integer.

7. the computational methods of the spiral CT image reconstruction mesoporous weighting that any one as claimed in claim 1 is accelerated based on GPU, it is characterised in that the aperture weighted value after described normalizationIt is stored in the shared video memory of GPU.