CN103366389A - CT (computed tomography) image reconstruction method - Google Patents

Abstract

The invention relates to a CT (computed tomography) image reconstruction method, which comprises the steps that detected signals of a detector are preprocessed, and projection data are obtained according to a detection result of the detector; the projection data are weighted according to scanned boundary parameters, and weighted projection data are subjected to one-dimensional convolutional filtering; the filtered projection data are mapped from a detector plane onto a virtual detector plane positioned in a rotation center according to scanned geometric parameters; filtered projection data, which are positioned on the virtual detector plane, of points to be reconstructed, are obtained through interpolation processing, and are subjected to three-dimensional back projection so as to obtain an initial image; and the initial image is used as an initial image of an expectation-maximization (EM) algorithm of ordered subsets, and is subjected to iterative reconstruction processing, and an imaged obtained after the iteration is terminated is output as a reconstructed image. By adopting the CT image reconstruction method, the CT image can be accurately reconstructed, and the definition of the reconstructed image can be improved.

Description

The CT image rebuilding method

Technical field

The present invention relates to image rebuilding method, especially a kind of CT(Computed Tomography) image rebuilding method.

Background technology

X ray computer fault imaging (CT, Computed Tomography) is interior inspection survey technology commonly used in the modern medicine, the CT technology mainly utilizes x-ray bundle and detector around human body rotating, and carry out continuously profile scanning, input computing machine after receiving the decay x-ray information of passing behind the human body by detector in each scanning process, carry out image reconstruction through robot calculator according to the decay x-ray information that receives, to obtain the image at human detection position.

Image reconstruction in the CT technology is a key factor that affects testing result, and is how clear, reappear to detect the image at position exactly, is the problem that those skilled in the art study always.

Summary of the invention

Provide hereinafter about brief overview of the present invention, in order to basic comprehension about some aspect of the present invention is provided.Should be appreciated that this general introduction is not about exhaustive general introduction of the present invention.It is not that intention is determined key of the present invention or pith, neither be intended to limit scope of the present invention.Its purpose only is that the form of simplifying provides some concept, with this as the in greater detail preorder of discussing after a while.

An object of the present invention is to provide a kind of can CT image rebuilding method clear, exactly reproduced image.

For achieving the above object, the invention provides a kind of CT image rebuilding method, comprising:

Detectable signal to detector carries out pre-service, and obtains data for projection according to the result of detection of detector;

Boundary parameter according to scanning is weighted described data for projection, and the data for projection after the weighting is carried out the one dimension convolutional filtering;

According to the scan geometry parameter, filtered data for projection is mapped to the dummy detector plane that is positioned at rotation center from the detector curved surface;

By interpolation processing obtain treat reconstruction point on the dummy detector plane through filtered data for projection, and carry out 3 D back projection, obtain initial pictures;

With the initial pictures of described initial pictures as the ordered subset expectation maximization value-based algorithm, carry out iterative approximation and process, the image after the output iteration stops is as reconstructed image.

The present invention by detectable signal is carried out pre-service, to data for projection be weighted, filtering etc., the error of the data for projection that recoverable obtains, can set up the higher initial pictures of sharpness and accuracy, in conjunction with order subset greatest hope value-based algorithm initial pictures being carried out iterative approximation processes, can accurately revise initial pictures, further promote sharpness and the accuracy of reconstructed image, for subsequent applications provides reliable basis.

Description of drawings

In order to be illustrated more clearly in the embodiment of the invention or technical scheme of the prior art, the below will do to introduce simply to the accompanying drawing of required use in embodiment or the description of the Prior Art, apparently, accompanying drawing in the following describes only is some embodiments of the present invention, for those of ordinary skills, under the prerequisite of not paying creative work, can also obtain according to these accompanying drawings other accompanying drawing.

The process flow diagram of the CT image rebuilding method that Fig. 1 provides for the embodiment of the invention;

Fig. 2 is the volume drawing physical model synoptic diagram that the present invention adopts when reconstructed image is carried out volume drawing.

Embodiment

For the purpose, technical scheme and the advantage that make the embodiment of the invention clearer, below in conjunction with the accompanying drawing in the embodiment of the invention, technical scheme in the embodiment of the invention is clearly and completely described, obviously, described embodiment is the present invention's part embodiment, rather than whole embodiment.The element of describing in an accompanying drawing of the present invention or a kind of embodiment and feature can combine with element and the feature shown in one or more other accompanying drawing or the embodiment.Should be noted that for purpose clearly, omitted expression and the description of parts that have nothing to do with the present invention, known to persons of ordinary skill in the art and processing in accompanying drawing and the explanation.Based on the embodiment among the present invention, the every other embodiment that those of ordinary skills obtain under the prerequisite of not paying creative work belongs to the scope of protection of the invention.

The process flow diagram of the image rebuilding method that Fig. 1 provides for the embodiment of the invention one.As shown in Figure 1, the method comprises:

Step S101: the detectable signal to detector carries out pre-service, and obtains data for projection according to the result of detection of detector; Can for example obtain the attenuation results of X ray behind material absorbing according to result of detection in this step and obtain data for projection, material was to its degree of absorption when this data for projection for example can be used for representing the X-ray scanning material.

Step S103: the boundary parameter according to scanning is weighted described data for projection, and the data for projection after the weighting is carried out the one dimension convolutional filtering; Described boundary parameter can comprise, for example reference position of frame and final position, the time shutter of bulb (ray tube), the useful signal unit scope of detection etc.In this step, can be weighted the data for projection p'(α that obtains after the weighting, beta, gamma to data for projection p (α, beta, gamma) by following formula (1)).

p'(α,β,γ)＝ω _2d(β,γ)·cosγ·p(α,β,γ) （1）

ω wherein _2dBe the one dimensional fourier transform function of object function 2d, α is wave vector α, and β is wave vector, and it is the point source of medium medium wave propagation for γ.

Can carry out the one dimension convolutional filtering to the data for projection after the weighting by following formula (2):

$p^{\′\′}(\mathrm{\α},\mathrm{\β},\mathrm{\γ})=\cos \mathrm{\α}\·p^{\′}(\mathrm{\α},\mathrm{\β},\mathrm{\γ})\⊗h\left(\mathrm{\γ}\right)---\left(2\right)$

Wherein, h (γ) is filter function.

Step S105: according to the scan geometry parameter, filtered data for projection is mapped to the dummy detector plane that is positioned at rotation center from the detector curved surface; This scan geometry parameter can comprise, the parameters such as the central passage during such as CT scan, radius of turn.

Step S107: by interpolation processing obtain treat reconstruction point on the dummy detector plane through filtered data for projection, and carry out 3 D back projection, obtain initial pictures;

This step is concrete according to treating that respectively reconstruction point carries out 3 D back projection through filtered data for projection on the dummy detector plane, obtains initial pictures.

If what work as pre-treatment is the data for projection of last scanning, then can obtain this initial pictures, otherwise, can repeat step S101, S103, S105, S107, until the data for projection of last scanning is also processed.

Step S109: with the initial pictures of described initial pictures as ordered subset expectation maximization value (OSEM) algorithm, carry out iterative approximation and process, the image after the output iteration stops is as reconstructed image.

Alternatively, step S109 comprises: 1), described data for projection is divided into n subset s ₁To s _n, repeat following steps, until after reconstructed image reached predefined termination criteria, the image when the output iteration stops was as described reconstructed image;

2), make

Wherein k represents the number of times of iteration, and its initial value is 0, when k=0,

Be described initial pictures, f ¹Be initial pictures when iterations k is 1;

3), by the theoretical data for projection of following formula (3) subset of computations:

${\mathrm{\μ}}_t^i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_{\mathrm{tj}}{f}_j^i,t\∈s_i---\left(3\right)$

Wherein,

J pixel of expression expression is to the image attenuation coefficient average of t bar projection ray, w _TjRepresent that j pixel is to the weighing factor of t bar projection ray, s _iRepresent the i subset, N represents number of pixels,

Represent that j pixel is to the initial pictures of t bar projection ray;

4), carrying out pixel by following formula (4) upgrades:

$f_j^{i+1}=f_j^i\frac{\underset{t\∈s_i}{\mathrm{\Σ}}\frac{p_t{w}_{\mathrm{tj}}}{{\mathrm{\μ}}_t^i}}{\underset{t\∈s_i}{\mathrm{\Σ}}{w}_{\mathrm{tj}}},(j=1,...,N)---\left(4\right)$

Wherein, p _tThe measured value that represents t bar projection ray;

5), with the theoretical data for projection of all subsets initial value as next iteration, namely Wherein, f ^N+1Represent the theoretical data for projection of all subsets as the initial value image of next iteration, namely all n subset projected datasets are the 1st the initial projected image of iteration on basis.

Such scheme by detectable signal is carried out pre-service, to data for projection be weighted, filtering etc., the error of the data for projection that recoverable obtains, can set up the higher initial pictures of sharpness and accuracy, in conjunction with order subset greatest hope value-based algorithm initial pictures being carried out iterative approximation processes, can accurately revise initial pictures, further promote sharpness and the accuracy of reconstructed image, for subsequent applications provides reliable basis.

Below in conjunction with this programme the ordered subset expectation maximization value-based algorithm is described in detail:

According to the physical geometry model that geometry and the scan protocols of CT frame are set up the CT system, the projection matrix W={w of computing system _Tj;

Initialization iterations k=0 sets non-negative initial pictures In this programme, adopt the initial pictures that obtains among the step S107;

Repeat the following step, until reconstructed image reaches predefined termination criteria, selectable standard comprises: iterations reaches the error of upper limit K, theoretical data for projection and actual measurement data for projection less than given threshold epsilon _p, the image update value

Less than given threshold epsilon _f

Enter the subset computing $f^1={\hat{f}}^k,k=k+1;$

From subset s ₁To s _n, at first carry out theoretical projection value and calculate, namely carry out the computing of formula (3).

Then carrying out pixel according to formula (4) upgrades.

The result that all subsets are calculated is as the initial value of next iteration, namely

The output iteration stops image

Reconstructed results as the OSEM algorithm.

Alternatively, in step S101, can comprise dark current removal, iterative filtering, light beam standardization, air calibration, negative logarithm operation, auto adapted filtering, water jet hardening correcting, the calibration of water mould etc.

The flow process that the present invention carries out the dark current removal is as follows:

Each probe unit to detector carries out signal sampling before exposure, with the value of sampled signal as dark current;

With when exposure to the sampled signal of each probe unit of detector with obtain the removal dark current after corresponding dark current subtracts each other after sampled signal.

For example, before the exposure, master controller sends dark current and removes instruction to the DAS(data acquisition system (DAS)), DAS starts dark current removal program after receiving instruction.To each probe unit, gather in a period of time the DAS signal of (such as 1s etc.), calculate its average as dark current, DAS collects during with exposure signal and dark current subtract each other, the signal after the removal of output dark current.

In the reality because the existence of dark current, DAS is not 0 at the signal that bulb does not have to collect in the situation of exposure yet, there is the stack of dark current in the signal that DAS collects when therefore exposing at every turn and produces error, the present invention by the way, can reduce the error of signals collecting, improve the accuracy of image reconstruction.

The present invention carries out iterative filtering by following formula (5) to the sampled signal of detector:

$x\left(\mathrm{k\Δt}\right)=\frac{y\left(\mathrm{k\Δt}\right)-\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\β}}_n{e}^{-\frac{\mathrm{\Δt}}{{\mathrm{\τ}}_n}}{S}_{\mathrm{nk}}}{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\β}}_n}---\left(5\right)$

Wherein

$S_{\mathrm{nk}}=x\left[(k-1)\mathrm{\Δt}\right]+e^{-\frac{\mathrm{\Δt}}{{\mathrm{\τ}}_n}}{S}_{n(k-1)}$

${\mathrm{\β}}_n={\mathrm{\α}}_n(1-e^{-\frac{\mathrm{\Δt}}{{\mathrm{\τ}}_n}})$

Wherein, Δ t is the time interval of signal sampling, and x (k Δ t) is the real integration sampling signal of described detector, and y (k Δ t) is the integration sampling signal of the described detector of measurement, N is the main attenuance component number that comprises in the detector family curve, τ _nThe damping time constant of n composition, τ _nTime constant, α _nThat time constant is τ _nThe time the relative intensity of attenuance component, k represents iterations, β _nRepresent the set of a plurality of projection angle actual measurement numbers, e is stochastic variable.

By carrying out the iterative filtering of detector sampled signal with upper type, the data error that recoverable is caused by the initial luminous and afterglow effect of detector.When detector employing rare-earth ceramic detector, it is initial luminous and impact afterglow effect is smaller, can consider to process with the attenuation model of low order.

The present invention carries out the light beam standardization and may further comprise the steps:

Arrange one at an end of detector and directly accept the X-ray irradiation with reference to probe unit, the intensity of real time record X-ray bundle is designated as I _R

With the realization standardization of being divided by of the signal of all the other probe units and signal with reference to probe unit.

Carry out by the way the light beam standardization, can remove the error that the undulatory property that produces owing to intensity time of the light that sends from bulb produces.

Based on above-mentioned light beam standardization, the present invention may further comprise the steps when carrying out air calibration:

Empty sampling step: in the situation of not putting into scanning object, carry out the bulb exposure and the detectable signal of detector is sampled;

The normalized gain obtaining step: with the sampled signal of each probe unit of described detector with reference to the sampled signal of the probe unit normalized gain divided by each probe unit that obtains detector;

Comparison step: in the Preset Time section, make detector from a plurality of angle sampled signals, obtain the mean value of a plurality of normalized gains of each probe unit, mean value and the predetermined threshold value of the normalized gain of each probe unit are compared, if the mean value of one of them normalized gain is then changed corresponding probe unit with not in default threshold range;

Look-up table forms step: the logarithm that calculates each gain mean value forms gain lookup;

Under the multiple condition of scanning, repeat described empty sampling step, normalized gain obtaining step, comparison step and look-up table and form step, form a plurality of gain lookup.

The frequency of air calibration for example, but is done 2 to 4 times according to place to use and temperature variation and decide every day.By air calibration, the detector gain of each passage of adjustable.

In a kind of example, can at first set correcting condition, comprise focus size, tube voltage, tube current and scanning bed thickness.Carry out the bulb exposure under the condition of not putting into scanning object, DAS begins image data simultaneously.

With the signal of each probe unit and signal with reference to the probe unit calculating normalized gain that is divided by.

The gain mean value of each probe unit (m, n) of calculated signals of a plurality of projection angles of repeated acquisition within a period of time (1s)

Itself and predefined detector normal operation gain threshold are compared, if gain mean value

Not in normal operation gain threshold scope, then this probe unit may lose efficacy needs to change.

The logarithm of calculated gains mean value Form gain lookup to accelerate follow-up negative logarithm operation.

For each available condition of work (tube voltage, scanning bed thickness, X-ray focus size), repeat above-mentioned steps, form a plurality of gain lookup.

The present invention specifically by the negative logarithm operation in the following formula (6), obtains data for projection:

$p(m,n)=-\ln \left(I_{\mathrm{out}}^N(m,n)\right)+\ln \left(I_{\mathrm{in}}(m,n)\right)$

Wherein, P (m, n) is data for projection, I _InThe light intensity of (m, n) emergent ray,

It is the light intensity of the ray that detector is accepted after the air calibration.

The present invention may further comprise the steps when carrying out auto adapted filtering:

Size according to following formula (7) calculation of filtered nuclear;

$D\left(x\right)\left\{\begin{array}{cc}\frac{2\mathrm{\β\λ}}{2\mathrm{\λ}+\mathrm{\β}(p(m,n)-T)}& \mathrm{if}(p(m,n)>T)\\ 0& \mathrm{otherwise}\end{array}\right.---\left(7\right)$

Wherein, D (x) is filtering core, β=15, and λ=75, T is the threshold value for the transmitted intensity of specific condition of scanning setting; T is an empirical parameter, and concrete value need to be determined by experiment.D (x) 〉=the 2nd, filtering core diameter choice function is along with drop shadow intensity increases and increases, if D (x)＜2 then data for projection is not carried out filtering.

If D (x) more than or equal to preset value, then carries out one-dimensional filtering according to following formula (8) to data for projection, particularly, data for projection p (m, n) is carried out along the detector one-dimensional filtering in D (x) the neighborhood scope radially:

$p_F(m,n)=\underset{|m_l-m|\≤D\left(x\right)}{\mathrm{\Σ}}p(m_l,n)\·\mathrm{\Λ}\left(\frac{m_l-m}{D\left(x\right)}\right)$

Wherein, p _F(m, n) is filtered data for projection, p (m _l, be the front data for projection of filtering n), Λ (x) is trigonometric function, m, n are the address corresponding to storer, m _lBe the arbitrary calculated value of interpolation.

For example, the X-ray line is followed Compound Poisson Distributions, and X-ray signal its noise ratio smaller for line is larger, may cause streak artifact.The present invention is in order to eliminate the streak artifact that is caused by noise, adopt adaptive filter method that data for projection is proofreaied and correct, the larger signal of X-ray line is minimized filtering or do not use filtering, and the signal smaller to the X-ray line then increases filtering, can reduce pseudo-shadow.

Before calculation of filtered nuclear size, can comprise focus size, tube voltage, tube current and scanning bed thickness etc. for each condition of scanning commonly used, corresponding X-ray beam intensity threshold value is set, for example, threshold value T.

The present invention specifically comprises when carrying out the water jet hardening correcting: carry out following steps under each tube voltage of ray tube:

11) calculate the length of passing the path of ray according to following formula (8):

l＝G ^-1(p)

G(l)＝∫μ(E)·l·f(E)dE （8）

Wherein, U is the tube voltage of ray tube, and l is the length of passing the path of ray tube emergent ray, and f (E) is the spectral distribution of ray tube emergent ray, and μ (E) is the attenuation coefficient of pure water, and p is the data for projection of actual measurement, and E is radiant force;

12) according under the given tube voltage U, the mapping relations of the length l of passing the path of ray tube emergent ray and the data for projection of actual measurement are set up look-up table (p ^-l) _U

13) obtain data for projection behind the water jet hardening correcting according to following formula (9):

p'＝l·μ(E) （9）

Wherein p ' is the data for projection behind the water jet hardening correcting.

Under each tube voltage of ray tube, carry out above-mentioned steps, can form a plurality of look-up tables, accelerate correction rate.

For example, the X-ray that sends from bulb has covered a spectral coverage that scope is very wide, and the projection value that causes recording no longer satisfies linear relationship with the path of passing material after bearing logarithm operation.The present invention considers that organ and glassware for drinking water that human body is main have the attenuation characteristic that relatively approaches, employing is based on the beam hardening correction of water, and according to above-mentioned correcting scheme, can obtain preferably calibration result, promote speed and the accuracy of image reconstruction, help to obtain clearly CT reconstructed image.

The present invention may further comprise the steps when carrying out the calibration of water mould:

21) under different sweep parameters, gather the data for projection of standard water mould, described standard water mould is placed on the center of scan vision, during concrete operations, can pass through the auxiliary of laser localized light comparing, the standard water mould is placed on the center of scan vision FOV.This sweep parameter can comprise the parameters such as tube current, tube voltage.

21) calculate μ ' according to following formula (9)

${\mathrm{\μ}}^{\′}={\mathrm{\μ}}_{\mathrm{air}}+(\mathrm{\μ}-\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{air}}})\·\frac{{\mathrm{\μ}}_{\mathrm{water}}-{\mathrm{\μ}}_{\mathrm{air}}}{\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{water}}}-\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{air}}}}---\left(9\right)$

Wherein,

For take the reconstructed image center as the center of circle, the water mould image attenuation coefficient average in the area-of-interest of radius 10cm,

Be the water mould image attenuation coefficient average in the air section, μ is for being the image attenuation coefficient average that does not have through the water jet hardening correcting, and μ ' is through the image attenuation coefficient average behind the water jet hardening correcting;

21) calculate the offset correction factor and convergent-divergent correction factor according to following formula (10):

${\mathrm{\μ}}_{\mathrm{scale}}=\frac{{\mathrm{\μ}}_{\mathrm{water}}-{\mathrm{\μ}}_{\mathrm{air}}}{\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{water}}}-\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{air}}}}$

${\mathrm{\μ}}_{\mathrm{offset}}=(\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{air}}}\·{\mathrm{\μ}}_{\mathrm{scale}}-{\mathrm{\μ}}_{\mathrm{air}})\·L---\left(10\right)$

Wherein, L is the distance from the ray tube focus to detector, μ _OffsetBe the offset correction factor, μ _ScaleBe the convergent-divergent correction factor;

21) obtain data for projection after water mould calibration according to following formula (11):

p″＝p′·μ _scale-μ _offset （11）

Wherein, " be the data for projection after the calibration of water mould, p ' is lower for the sweep parameter in correspondence arranges, and individual data for projection by gathering also can be the data for projection behind the water jet hardening correcting for p.

The CT value of reconstructed image Air should be-1000HU in theory, the CT value of pure water should be 0HU, the possibility of result of actual reconstruction has certain deviation, the present invention adopts the standard water mould, by with upper type data for projection being proofreaied and correct, can realize the CT value deviation of reconstructed image is revised, improve the accuracy of image reconstruction.Alternatively, determine emending frequency according to the bulb time shutter, the calibration of water mould can general three months be done once to half a year, upgrades correction parameter.

Alternatively, in step S109, obtain reconstructed image after, also can comprise the step of reconstructed image being carried out aftertreatment, for example, bone beam hardening correction, the removal of pseudo-shadow, visualization processing etc., the below is introduced it respectively.

The present invention can carry out the bone beam hardening correction behind the output reconstructed image, specifically comprise:

Under the tube voltage of predefined ray tube, according to the spectral distribution of emergent ray, the attenuation coefficient of bone and the data for projection that passes the path calculating bone of bone;

Data for projection to bone carries out the water jet hardening correcting, the data for projection behind the acquisition water jet hardening correcting;

The data for projection of the method that adopts additional second order correction after with the water jet hardening correcting is mapped to desirable bone response curve, record second order correction parameter;

Image after rebuilding is carried out Threshold segmentation, isolate bone tissue;

Bone tissue is simulated forward projection, to produce the data for projection that only contains bone;

The data for projection that only contains bone is squared, generate the projection error data;

According to the projection error data, reconstruct an image that only contains pseudo-shadow;

The image that only contains pseudo-shadow be multiply by the second order correction parameter, from reconstructed image, deduct the image that only contains pseudo-shadow, the image after the beam alignment of generation bone.

By carrying out the bone beam hardening correction with upper type, can reduce because the pseudo-shadow that the marked difference of attenuation characteristic causes between bone tissue and the water tissue sharpness of lifting reconstructed image.

The present invention can carry out pseudo-shadow Transformatin behind the output reconstructed image, specifically comprise:

31) reconstructed image is transformed to polar coordinate system (r, θ) by theorem in Euclid space coordinate system (x, y), the ring in the theorem in Euclid space is mapped as the straight line under the polar coordinates at this moment;

32) under polar coordinate system, according to following formula (12) each pixel is carried out radially one-dimensional filtering: particularly, the filtering core k (r) that adopts 1 * N carries out one-dimensional filtering on the r direction to each pixel p (r, θ);

$f(r,\mathrm{\θ})=|k\left(r\right)\⊗p(r,\mathrm{\θ})|---\left(12\right)$

Wherein, p (r, θ) is pixel value, and k (r) is filtering core, and f (r, θ) is filtered pixel value;

33) if Normal＜f (r, θ)＜Edge, wherein Normal is the convolution threshold value of normal pixel, Edge is the convolution threshold value of image normal boundary, then pixel value p (r, θ) is judged as ring artifact; The value of Normal and Edge need to be by testing to obtain empirical parameter setting to different tissues, different reconstruction filter kernel;

34) according to following formula (13), adopt the pixel value in the radial neighborhood that the pixel value that is judged to be ring artifact is carried out match:

$p(r,\mathrm{\θ})=\underset{r_i\∈[r-\mathrm{\Δr},r+\mathrm{\Δr}]}{\mathrm{\Σ}}{\mathrm{\α}}_i\·p(r_i,\mathrm{\θ})---\left(13\right)$

α wherein _iBe wave vector;

35) the theorem in Euclid space coordinate system is changed in the polar coordinate system image inversion after the match, eliminated ring artifact;

36) strong attenuation occurs when X-ray penetrating metal object, the X-ray subnumber that detector receives is deficient, thereby signal noise causes occurring streak artifact in the image greatly, also can carry out following steps in this step and eliminate strip artifact:

Reconstructed image is adopted Threshold segmentation, be partitioned into the metal object zone;

37) its corresponding original projection data area [γ is calculated in each metal object zone _a, γ _b], carry out following correction to obtain revised data for projection p for the data for projection p (α, beta, gamma) of each projection angle β _F(α, beta, gamma):

$p_F(\mathrm{\α},\mathrm{\β},\mathrm{\γ})=\left\{\begin{array}{cc}p(\mathrm{\α},\mathrm{\β},\mathrm{\γ})& \mathrm{\γ}\∉[{\mathrm{\γ}}_a,{\mathrm{\γ}}_b]\\ \frac{{\mathrm{\γ}}_b-\mathrm{\γ}}{{\mathrm{\γ}}_b-{\mathrm{\γ}}_a}p(\mathrm{\α},\mathrm{\β},{\mathrm{\γ}}_a)+\frac{\mathrm{\γ}-{\mathrm{\γ}}_a}{{\mathrm{\γ}}_b-{\mathrm{\γ}}_a}p(\mathrm{\α},\mathrm{\β},{\mathrm{\γ}}_b)& \mathrm{\γ}\∈[{\mathrm{\γ}}_a,{\mathrm{\γ}}_b]\end{array}\right.---\left(14\right)$

38) adopt revised data for projection to carry out image reconstruction, obtain the not image of containing metal object information;

A gray scale drawing coefficient is multiply by in the metal object zone that is partitioned into, merge mutually the image that obtains after metal artifacts is proofreaied and correct with the image of containing metal object information not; Gray scale drawing coefficient s ₂Be an empirical parameter, need to be determined by experiment.

Alternatively, behind the output reconstructed image, also can carry out visualization processing to this reconstructed image, specifically describe as follows:

To the CT image after rebuilding, can carry out mutual visual of two and three dimensions by computing machine, with help the doctor more directly perceived, grasp patient's physiology, pathological information all sidedly.The present invention program intends providing following concentrated visualization technique: in real time section is browsed, many planar reconstruction and volume drawing.

In real time section is browsed

An in real time section is set in rebuilding main interface browses window, show in real time reconstruction progress and current reconstruction section, and the selection of rebuilding section function of browse is provided.

Many planar reconstruction

Set up one based on patient's coordinate system, three orthogonal planes are arranged: coronal-plane, sagittal plane, axle shape face, any and uneven plane, these planes is called as the inclined-plane.By in different directions to the volumetric image interpolation, produce reconstructed image.Drag the plane of choosing and to realize that many planes of three dimensional CT volume data show and browse.

The key step that many planar reconstruction show is as follows:

Set up many planar reconstruction displaing coordinate system, comprising: world coordinate system, object coordinates system, eye coordinate, device coordinate system etc.;

For the each time operation of image rotation, translation and convergent-divergent, the Renewal model transformation matrix is mapped to world coordinate system with summit, 3-D display zone from object coordinates system, and the three-dimensional CT images data are tied on the viewing area in the world coordinates as texture;

By the plane that user interactions select to determine need shows, calculate each position in world coordinate system, the summit that needs display plane, then obtain 2 d texture on the display plane at the enterprising row interpolation of three-D grain;

With 2 d texture and display plane binding, upgrade the view transformation matrix with display plane from the world coordinate system projection mapping to eye coordinate, then carry out device transform display plane is mapped to device coordinate system, in drawing window, carry out at last many planar reconstruction and show.

Volume drawing

As shown in Figure 2, the Volume Rendering Techniques of the present invention's proposition at first needs to select an observer visual angle and an imaginary screen.Then one group of virtual ray penetrates object along view directions, ray can be dispersed or be parallel, by a mapping function, the intensity of a pixel or color are determined by predefined mapping relations by all pixels that intersect with ray in the three-dimensional body on the imaginary screen.Intending adopting following steps to carry out volume drawing in the scheme shows:

Set up a volume drawing displaing coordinate system, comprising: world coordinate system, object coordinates system, eye coordinate, device coordinate system etc.;

Set volume drawing and show the transport function that adopts, describe the mapping relations between the information such as pixel value, pixel gradient and the corresponding optical properties (opacity, color etc.);

For the each time operation of image rotation, translation and convergent-divergent, the Renewal model transformation matrix is mapped to world coordinate system with the three-dimensional CT images data from object coordinates system;

Launch a light from each pixel of projection plane, pass the three-dimensional CT images data, calculate the light attribute of each sampled point S.For each sampled point S, not only self have certain density and color, also can be with light reflection on every side to viewpoint, the two has consisted of the emergent ray C (s) that S is ordered, and can be expressed as:

C(s)＝Ka*Ca+Kd*Cl*Co(s)(N(s)·L(s))+Ks*Cl*(N(s)·H(s)) ^ns (15)

Wherein Ca is surround lighting, Ka is the surround lighting coefficient, Kd is diffuse-reflection factor, Cl is light source colour, and Co (s) is the color of S point itself, the gradient vector that N (s) is ordered for S, L (s) is light source vector, Ks is specularity factor, and H (s) is high light vector, and ns is high optical index.

Because emergent ray also can be subject to the barrier effect of other point the process that projects viewpoint from S point, so the S point is the pad value of C (s) to the light contribution of viewpoint.Because each particle the Ray(ray) can have a light contribution margin to viewpoint, therefore suppose that the light total amount that viewpoint receives in the Ray direction is I, then I can be expressed as the form of a line integral:

$I={\∫}_0^{R}C\left(s\right)*\mathrm{\μ}\left(s\right)*\exp (-{\∫}_0^s\mathrm{\μ}\left(t\right)\mathrm{dt})\mathrm{ds}---\left(16\right)$

Wherein μ (s) is the density value of sampled point S, and R is the length of Ray.

Order according to Front-to-back(from front to back) is carried out the mixing of sampled point light attribute, suppose from viewpoint, on Ray every point of Δ S distance samples, then formula (16) can approximately equivalent in the mixed formulation of a recurrence:

$\left\{\begin{array}{c}c=C(i\·\mathrm{\ΔS})*\mathrm{\α}(i\·\mathrm{\ΔS})*(1-\mathrm{\α})+c\\ \mathrm{\α}=\mathrm{\α}(i\·\mathrm{\ΔS})*(1-\mathrm{\α})+\mathrm{\α}\end{array}\right.$

Wherein c is the light color accumulated value of viewpoint, and α is the opacity accumulated value of viewpoint, and C (i Δ S) is the emergent ray of i sample point, and α (i Δ S) is the opacity of i sample point.

Generate the final two-dimensional projection image that shows, the color value of corresponding each pixel is that c is to the ratio of α.

In the various embodiments described above of the present invention, the sequence number of embodiment only is convenient to describe, and does not represent the quality of embodiment.Description to each embodiment all emphasizes particularly on different fields, and does not have the part of detailed description among certain embodiment, can be referring to the associated description of other embodiment.

One of ordinary skill in the art will appreciate that: all or part of step that realizes said method embodiment can be finished by the relevant hardware of programmed instruction, aforesaid program can be stored in the computer read/write memory medium, this program is carried out the step that comprises said method embodiment when carrying out; And aforesaid storage medium comprises: the various media that can be program code stored such as ROM (read-only memory) (Read-Only Memory is called for short ROM), random access memory (Random Access Memory is called for short RAM), magnetic disc or CD.

In the embodiment such as apparatus and method of the present invention, obviously, after can decomposing, make up and/or decompose, each parts or each step reconfigure.These decomposition and/or reconfigure and to be considered as equivalents of the present invention.Simultaneously, in the above in the description to the specific embodiment of the invention, can in one or more other embodiment, use in identical or similar mode for the feature that a kind of embodiment is described and/or illustrated, combined with the feature in other embodiment, or the feature in alternative other embodiment.

Should emphasize that term " comprises/comprise " existence that refers to feature, key element, step or assembly when this paper uses, but not get rid of the existence of one or more further feature, key element, step or assembly or additional.

It should be noted that at last: although above the present invention and the advantage thereof of having described in detail is to be understood that and can carries out various changes, alternative and conversion in the situation that does not exceed the spirit and scope of the present invention that limited by appended claim.And scope of the present invention is not limited only to the specific embodiment of the described process of instructions, equipment, means, method and step.The one of ordinary skilled in the art will readily appreciate that from disclosure of the present invention, can use according to the present invention carry out with the essentially identical function of corresponding embodiment described herein or obtain result essentially identical with it, existing and want exploited process, equipment, means, method or step future.Therefore, appended claim is intended to comprise such process, equipment, means, method or step in their scope.

Claims (10)

1. a CT image rebuilding method is characterized in that, comprising:

Detectable signal to detector carries out pre-service, and obtains data for projection according to the result of detection of detector;

Boundary parameter according to scanning is weighted described data for projection, and the data for projection after the weighting is carried out the one dimension convolutional filtering;

According to the scan geometry parameter, filtered data for projection is mapped to the dummy detector plane that is positioned at rotation center from the detector curved surface;

By interpolation processing obtain treat reconstruction point on the dummy detector plane through filtered data for projection, and carry out 3 D back projection, obtain initial pictures;

With the initial pictures of described initial pictures as the ordered subset expectation maximization value-based algorithm, carry out iterative approximation and process, the image after the output iteration stops is as reconstructed image.

2. CT image rebuilding method according to claim 1 is characterized in that, the detectable signal of detector is carried out pre-service specifically comprise:

Each probe unit to detector carries out signal sampling before exposure, with the value of sampled signal as dark current;

With when exposure to the sampled signal of each probe unit of detector with obtain the removal dark current after corresponding dark current subtracts each other after sampled signal.

3. CT image rebuilding method according to claim 1 is characterized in that, the detectable signal of detector is carried out pre-service specifically comprise:

By following formula the sampled signal of described detector is carried out iterative filtering:

$x\left(\mathrm{k\Δt}\right)=\frac{y\left(\mathrm{k\Δt}\right)-\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\β}}_n{e}^{-\frac{\mathrm{\Δt}}{{\mathrm{\τ}}_n}}{S}_{\mathrm{nk}}}{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\β}}_n}$

Wherein

$S_{\mathrm{nk}}=x\left[(k-1)\mathrm{\Δt}\right]+e^{-\frac{\mathrm{\Δt}}{{\mathrm{\τ}}_n}}{S}_{n(k-1)}$

${\mathrm{\β}}_n={\mathrm{\α}}_n(1-e^{-\frac{\mathrm{\Δt}}{{\mathrm{\τ}}_n}})$

Δ t is the time interval of signal sampling, and x (k Δ t) is the real integration sampling signal of described detector, and y (k Δ t) is the integration sampling signal of the described detector of measurement, and N is the main attenuance component number that comprises in the detector family curve, τ _nThe damping time constant of n composition, τ _nTime constant, α _nThat time constant is τ _nThe time the relative intensity of attenuance component, k represents iterations, β _nRepresent the set of a plurality of projection angle actual measurement numbers, e is stochastic variable.

4. CT image rebuilding method according to claim 1 is characterized in that, the detectable signal of detector is carried out pre-service specifically comprise:

Empty sampling step: in the situation of not putting into scanning object, carry out the ray tube exposure and the detectable signal of detector is sampled;

The normalized gain obtaining step: with the sampled signal of each probe unit of described detector with reference to the sampled signal of the probe unit normalized gain divided by each probe unit that obtains detector;

Comparison step: in the Preset Time section, make detector from a plurality of angle sampled signals, obtain the mean value of a plurality of normalized gains of each probe unit, mean value and the predetermined threshold value of the normalized gain of each probe unit are compared, if the mean value of one of them normalized gain is then changed corresponding probe unit with not in default threshold range;

Look-up table forms step: the logarithm that calculates each gain mean value forms gain lookup;

Under the multiple condition of scanning, repeat described empty sampling step, normalized gain obtaining step, comparison step and look-up table and form step, form a plurality of gain lookup.

5. CT image rebuilding method according to claim 4 is characterized in that, the detectable signal of detector is carried out pre-service specifically comprise:

Obtain data for projection according to following formula:

$p(m,n)=-\ln \left(I_{\mathrm{out}}^N(m,n)\right)+\ln \left(I_{\mathrm{in}}(m,n)\right)$

Wherein, p (m, n) is data for projection, I _InThe light intensity of (m, n) emergent ray,

It is the light intensity of the ray that detector is accepted after the air calibration.

6. CT image rebuilding method according to claim 5 is characterized in that, the detectable signal of detector is carried out pre-service specifically comprise:

Size according to following formula calculation of filtered nuclear;

$D\left(x\right)\left\{\begin{array}{cc}\frac{2\mathrm{\β\λ}}{2\mathrm{\λ}+\mathrm{\β}(p(m,n)-T)}& \mathrm{if}(p(m,n)>T)\\ 0& \mathrm{otherwise}\end{array}\right.$

Wherein, D (x) is filtering core, β=15, and λ=75, T is the threshold value for the transmitted intensity of specific condition of scanning setting;

If D (x) then carries out one-dimensional filtering according to following formula to described data for projection more than or equal to preset value $p_F(m,n)=\underset{|m_l-m|\≤D\left(x\right)}{\mathrm{\Σ}}p(m_l,n)\·\mathrm{\Λ}\left(\frac{m_l-m}{D\left(x\right)}\right)$

Wherein, p _F(m, n) is filtered data for projection, p (m _l, be the front data for projection of filtering n), m, n are the address corresponding to storer, m _lBe the arbitrary calculated value of interpolation.

7. CT image rebuilding method according to claim 6 is characterized in that,

The detectable signal of detector is carried out pre-service specifically to be comprised: carry out following steps under each tube voltage of ray tube:

Calculate the length of passing the path of ray according to following formula:

l＝G ^-1(p)

G(l)＝∫μ(E)·l·f(E)dE

Wherein, U is the tube voltage of ray tube, and l is the length of passing the path of ray tube emergent ray, and f (E) is the spectral distribution of ray tube emergent ray, and μ (E) is the attenuation coefficient of pure water, and p is the data for projection of actual measurement, and E is radiant force;

Under given tube voltage U, the mapping relations of the length l of passing the path of ray tube emergent ray and the data for projection of actual measurement are set up look-up table (p-l) _U

Obtain data for projection behind the water jet hardening correcting according to following formula:

P'=l μ (E), wherein p' is the data for projection behind the water jet hardening correcting;

And/or,

The detectable signal of detector is carried out pre-service specifically to be comprised:

Under different sweep parameters, gather the data for projection of standard water mould, described standard water mould is placed on the center of scan vision;

Calculate μ ' according to following formula

${\mathrm{\μ}}^{\′}={\mathrm{\μ}}_{\mathrm{air}}+(\mathrm{\μ}-\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{air}}})\·\frac{{\mathrm{\μ}}_{\mathrm{water}}-{\mathrm{\μ}}_{\mathrm{air}}}{\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{water}}}-\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{air}}}}$

Wherein,

For take the reconstructed image center as the center of circle, the water mould image attenuation coefficient average in the area-of-interest of radius 10cm, Be the water mould image attenuation coefficient average in the air section, μ is not for there being the image attenuation coefficient average through the water jet hardening correcting, and μ ' is through the image attenuation coefficient average behind the water jet hardening correcting;

Calculate the offset correction factor and convergent-divergent correction factor according to following formula:

${\mathrm{\μ}}_{\mathrm{scale}}=\frac{{\mathrm{\μ}}_{\mathrm{water}}-{\mathrm{\μ}}_{\mathrm{air}}}{\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{water}}}-\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{air}}}}$

${\mathrm{\μ}}_{\mathrm{offset}}=(\stackrel{\‾}{{\mathrm{\μ}}_{\mathrm{air}}}\·{\mathrm{\μ}}_{\mathrm{scale}}-{\mathrm{\μ}}_{\mathrm{air}})\·L$

Wherein, L is the distance from the ray tube focus to detector, μ _OffsetBe the offset correction factor, μ _ScaleBe the convergent-divergent correction factor;

Obtain data for projection after water mould calibration according to following formula:

p″=p′·μ _scale-μ _offset

Wherein, p " is the data for projection after the calibration of water mould.

8. CT image rebuilding method according to claim 1 is characterized in that, also comprises after the output reconstructed image:

Under the tube voltage of predefined ray tube, according to the spectral distribution of emergent ray, the attenuation coefficient of bone and the data for projection that passes the path calculating bone of bone;

Data for projection to bone carries out the water jet hardening correcting, the data for projection behind the acquisition water jet hardening correcting;

The data for projection of the method that adopts additional second order correction after with the water jet hardening correcting is mapped to desirable bone response curve, record second order correction parameter;

Image after rebuilding is carried out Threshold segmentation, isolate bone tissue;

Bone tissue is simulated forward projection, to produce the data for projection that only contains bone;

The data for projection that only contains bone is squared, generate the projection error data;

According to the projection error data, reconstruct an image that only contains pseudo-shadow;

The image that only contains pseudo-shadow be multiply by the second order correction parameter, from reconstructed image, deduct the image that only contains pseudo-shadow, the image after the beam alignment of generation bone.

9. CT image rebuilding method according to claim 1 is characterized in that, also comprises after the output reconstructed image:

Reconstructed image is transformed to polar coordinate system (r, θ) by theorem in Euclid space coordinate system (x, y);

Under polar coordinate system, according to following formula each pixel is carried out radially one-dimensional filtering:

$f(r,\mathrm{\θ})=|k\left(r\right)\⊗p(r,\mathrm{\θ})|$

Wherein, p (r, θ) is pixel value, and k (r) is filtering core, and f (r, θ) is filtered pixel value

If Normal＜f (r, θ)＜Edge, wherein Normal is the convolution threshold value of normal pixel, and Edge is the convolution threshold value of image normal boundary, and then pixel value p (r, θ) is judged as ring artifact;

Adopt the pixel value in the radial neighborhood that the pixel value that is judged to be ring artifact is carried out match:

$p(r,\mathrm{\θ})=\underset{r_i\∈[r-\mathrm{\Δr},r+\mathrm{\Δr}]}{\mathrm{\Σ}}{\mathrm{\α}}_i\·p(r_i,\mathrm{\θ})$

The theorem in Euclid space coordinate system is changed in polar coordinate system image inversion after the match, eliminated ring artifact;

Reconstructed image is adopted Threshold segmentation, be partitioned into the metal object zone;

To whenever _OneIts corresponding original projection data area [γ is calculated in individual metal object zone _a, γ _b], carry out following correction to obtain revised data for projection p for the data for projection p (α, beta, gamma) of each projection angle β _F(α, beta, gamma):

$p_F(\mathrm{\α},\mathrm{\β},\mathrm{\γ})=\left\{\begin{array}{cc}p(\mathrm{\α},\mathrm{\β},\mathrm{\γ})& \mathrm{\γ}\∉[{\mathrm{\γ}}_a,{\mathrm{\γ}}_b]\\ \frac{{\mathrm{\γ}}_b-\mathrm{\γ}}{{\mathrm{\γ}}_b-{\mathrm{\γ}}_a}p(\mathrm{\α},\mathrm{\β},{\mathrm{\γ}}_a)+\frac{\mathrm{\γ}-{\mathrm{\γ}}_a}{{\mathrm{\γ}}_b-{\mathrm{\γ}}_a}p(\mathrm{\α},\mathrm{\β},{\mathrm{\γ}}_b)& \mathrm{\γ}\∈[{\mathrm{\γ}}_a,{\mathrm{\γ}}_b]\end{array}\right.$

Adopt revised data for projection to carry out image reconstruction, obtain the not image of containing metal object information;

A gray scale drawing coefficient is multiply by in the metal object zone that is partitioned into, merge mutually the image that obtains after metal artifacts is proofreaied and correct with the image of containing metal object information not.

10. CT image rebuilding method according to claim 1 is characterized in that, described initial pictures is carried out iterative approximation process, and specifically comprises:

Described data for projection is divided into n subset s ₁To s _n, repeat following steps, until after reconstructed image reached predefined termination criteria, the image when the output iteration stops was as described reconstructed image;

Make

K=k+1, wherein k represents the number of times of iteration, its initial value is 0, when k=0,

Be described initial pictures, f ¹Be initial pictures when iterations k is 1;

By the theoretical data for projection of following formula subset of computations:

${\mathrm{\μ}}_t^i=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{w}_{\mathrm{tj}}{f}_j^i,t\∈s_i$

Wherein,

Represent that j pixel is to the image attenuation coefficient average of t bar projection ray, w _TjRepresent that j pixel is to the weighing factor of t bar projection ray, s _iRepresent the i subset, N represents number of pixels,

Represent the i subset, the projected image of j pixel;

Carrying out pixel by following formula upgrades:

$f_j^{i+1}=f_j^i\frac{\underset{t\∈s_i}{\mathrm{\Σ}}\frac{p_t{w}_{\mathrm{tj}}}{{\mathrm{\μ}}_t^i}}{\underset{t\∈s_i}{\mathrm{\Σ}}{w}_{\mathrm{tj}}},(j=1,...,N)$

Wherein, p _tThe measured value that represents t bar projection ray;

With the theoretical data for projection of all subsets initial value as next iteration, namely

Wherein, f ^N+1Represent all n subset projected datasets and be the 1st the initial projected image of iteration on basis.

Images