CN105321155A - Ring artifact elimination method for CBCT image - Google Patents

Abstract

The invention relates to a ring artifact elimination method for a CBCT image and belongs to the field of medical image processing. According to the invention, firstly, the CBCT image with a ring artifact is transformed from a Cartesian coordinate system to a polar coordinate system; then the CT image P in the polar coordinate system is subjected to smoothing processing by adopting an L0 norm filtering method; next, a smooth image P' obtained after filtering is subtracted from an original image P in the polar coordinate system to obtain an artifact information image Pk with partial original detail information, then the Pk is subjected to extraction and removal of an artifact template by using an OSTU multi-threshold segmentation method and a median method, and the rest of detail information is compensated into the filtered image P'; and finally, the image P' compensated with the detail information is transformed to the Cartesian coordinate system, so that a cone-beam CT image without the ring artifact can be obtained. Compared with the prior art, when detail and edge information in the cone-beam CT image is reserved to the greatest extent, the ring artifact on the CT image, which influences image quality and reality, are effectively removed without obvious residues.

Description

A kind of CBCT image annular artifact eliminating method

Technical field

The present invention relates to a kind of elimination CBCT image annular artifact method, belong to field of medical image processing.

Background technology

In recent years, Conebeam computed tomography imaging system (Cone-beamComputedTomography) is widely used at medical fields such as clinical diagnosis, 3D implants.As the one of fault imaging, conical beam CT has the many merits such as imaging is accurate, correlation line use ratio is high.But, due to the restriction of the multiple factors such as system process, design and reconstruction algorithm, often have in conical beam CT image and be a series ofly the center of circle with reconstructed center and gray-scale value is different from the donut of surrounding pixel, i.e. annular artifact.The quality that have impact on image that the appearance of annular artifact is serious and authenticity, disturb further clinical diagnosis and process.Thus, how when not affecting image resolution ratio, the emphasis becoming current industry research to be removed to annular artifact.

Occurring so far from conical beam CT imaging system, there is the multiple antidote for annular artifact and has achieved certain effect in industry.In recent years, the removal work of annular artifact is mainly carried out and is roughly divided into two large classes on image, respectively:

(1) based on the pre-treating method of projection sinogram;

(2) based on the post-processing approach rebuilding image;

For first kind antidote, its cardinal principle is in the projection sinogram that obtains before image is rebuilt, and annular artifact is usually expressed as the straight line that on vertical direction, many are parallel to each other.Compared to donut, vertical striped has and is more easy to detected and the features such as elimination.2009, the people such as BeatMunch propose method wavelet decomposition and Fourier's low-pass filter combined, in order to the pseudo-shadow information and raw information that make image obtain stricter differentiation, first multi-scale wavelet decomposition is done to the projection sinogram of CT image, thus the vertical frequency band containing image in the vertical direction partial information under obtaining different scale.Thereafter, the Fast Fourier Transform (FFT) in spatial domain done to all vertical frequency bands and carries out low-pass filtering, just can obtain detailed information finally by the method for anti-wavelet decomposition and retain image after more intact process.Although this method to eliminate in low-pass filtering process filter function to a great extent for the impact of raw information, but as the common fault of pre-treating method, storage space shared by projection sinogram is huge, and the operation of this series of complex will expend the regular hour thus; Simultaneously for different pictures, also need to carry out artificial selection to multiple different parameter thus result in the decline of rebuilding efficiency in many aspects.Further, in some cases, this method may cause the generation of the residual of artifact and new artifact.

In order to save storage space, improving the efficiency that artifact is removed, the correction result of image being had and embodies more intuitively.In recent years, industry has engendered the multiple post-processing approach for rebuilding image.Similar to the pre-treating method based on projection sinogram, in order to change the geometrical property of annular artifact, reducing artifact and detecting and the difficulty in removal process, first need carry out polar coordinate transform to image after reconstruction.According to the mathematical property of this type of conversion, the ring artifact in cartesian coordinate system will show as stripe artifact in polar coordinate system, through follow-up a series of process and inverse coordinate transform, just can obtain the CT image reconstruction after annular artifact correction.The existing post-processing approach for rebuilding image have employed the filtering of different modes mostly at present, causes edge fog and the loss of detail of image thus.2007, the method based on unidirectional adaptive smooth process that the people such as Xia Xiongjun propose, eliminating in CT image the original detailed information as far as possible maintained while annular artifact in image, what reduce its edge is fuzzy, is a kind of excellent in the post-processing approach rebuilding image.But due to parameter selection with need consider that artifact removes the problems such as the balance retained with details, when annular artifact intensity transformation is larger, this kind of method can cause artifact to a certain extent residual.

Summary of the invention

The object of the invention is for problem existing in the post-processing approach at present with stronger development prospect, propose a kind of CBCT image annular artifact eliminating method, the method on the basis of post-processing approach, by L ₀norm filtering, OSTU Threshold segmentation combine with artifact template extraction, when retaining image border and detailed information, effectively remove the annular artifact in CT image.

The inventive method is achieved through the following technical solutions:

A kind of CBCT image annular artifact eliminating method, its protocol step is as follows:

Step one, transforms to artifact by original image I by cartesian coordinate system and shows as the image P that the polar coordinate system of parallel lines on vertical direction obtains under polar coordinate system;

As preferably, after obtaining image P described in step one, the most contiguous difference approach is adopted to compensate to reach resolution requirement to it.

Step 2, adopts L ₀the method of norm filtering obtains smoothed image P' to the smoothing process of image P under polar coordinate system;

Step 3, by subtracting each other the artifact hum pattern P obtaining having the original detailed information of part by the smoothed image P' obtained after the original image P under polar coordinate system and filtering _k, then to P _kextraction and the removal of carrying out artifact template obtain P _k', finally will eliminate artifact template and only comprise the P of detailed information _k' compensate back filtered image P ';

As preferably, to P _kthe extraction carrying out artifact template is completed by following process with removal:

(1) OSTU multi-threshold segmentation method is adopted to extract the area-of-interest having obvious stripe-shaped artifact;

As preferably, the described area-of-interest to having obvious stripe-shaped artifact is extracted and is completed by following process: to former figure P and artifact hum pattern P _kdo multi-threshold segmentation respectively and in area-of-interest, get its common factor and obtain final area-of-interest.

(2) at artifact hum pattern P _karea-of-interest in utilize median method to carry out extracting and remove to artifact template;

As preferably, the described median method that utilizes carries out extracting to artifact template and removes further by getting front-seat higher with the lower point of some gray-scale values that removes of intermediate value to obtain better effect.

Step 4, is converted into cartesian coordinate system by the image P ' that compensate for detailed information, just can obtain the conical beam CT image of a width ringless-type artifact.

Beneficial effect

The present invention applies in real tissue conical beam CT image, contrast with existing method, to retain in conical beam CT image while details and marginal information to the full extent, effectively eliminate the annular artifact that it affects picture quality and authenticity and there is no obvious remaining.Meanwhile, except the overall gray value of image there occurs except slight conversion, there is not obvious decline in the resolution of image itself, embodies superiority of the present invention yet.

Accompanying drawing explanation

Fig. 1 is that the one that proposes of the present invention is based on L ₀the basic procedure of the annular artifact minimizing technology of norm filtering.

Fig. 2 is the reconstruction image carried out before and after polar coordinate transform, and wherein (a) is brain CBCT original image F; B () is the brain image P under polar coordinate system.

The result design sketch of Fig. 3 to be the present invention with the comparatively classical annular artifact minimizing technology (wavelet-Fourierfiltering) based on wavelet decomposition and Fourier transform and the method (adaptivelowpassfilter) based on unidirectional adaptive smooth process carry out a width head conical beam CT image 1 annular artifact removal; Wherein, (a) is protocephalic region cone beam ct reconstruction image 1; B () is the image after the method process that proposes of the present invention; C () is the Yuan Tu local magnification region chosen; D () is the former figure of local magnification region; E () is the local area image after the process of adaptivelowpassfilter method; F () is the local area image after the process of wavelet-Fourierfiltering method; G () is the L that the present invention adopts ₀the result of smoothing method when Selecting parameter λ is 0.03.

The result design sketch of Fig. 4 to be the present invention with the comparatively classical annular artifact minimizing technology (wavelet-Fourierfiltering) based on wavelet decomposition and Fourier transform and the method (adaptivelowpassfilter) based on unidirectional adaptive smooth process carry out a width head conical beam CT image 2 annular artifact removal; Wherein, (a) is protocephalic region cone beam ct reconstruction image 2; B () is the image after the method process that proposes of the present invention; C () is the Yuan Tu local magnification region chosen; D () is the former figure of local magnification region; E () is the local area image after the process of adaptivelowpassfilter method; F () is the local area image after the process of wavelet-Fourierfiltering method; G () is the L that the present invention adopts ₀the result of smoothing method when Selecting parameter λ is 0.03.

The result design sketch of Fig. 5 to be the present invention with the comparatively classical annular artifact minimizing technology (wavelet-Fourierfiltering) based on wavelet decomposition and Fourier transform and the method (adaptivelowpassfilter) based on unidirectional adaptive smooth process carry out width neck conical beam CT image annular artifact removal; Wherein, (a) is former neck conical beam CT image; B () is the image after the method process that proposes of the present invention; C () is the Yuan Tu local magnification region chosen; D () is the former figure of local magnification region; E () is the local area image after the process of adaptivelowpassfilter method; F () is the local area image after the process of wavelet-Fourierfiltering method; G () is the L that the present invention adopts ₀the result of smoothing method when Selecting parameter λ is 0.005.

The result design sketch of Fig. 6 to be the present invention with the comparatively classical annular artifact minimizing technology (wavelet-Fourierfiltering) based on wavelet decomposition and Fourier transform and the method (adaptivelowpassfilter) based on unidirectional adaptive smooth process carry out width basis cranii conical beam CT image annular artifact removal; Wherein, (a) is cone beam ct reconstruction image at the bottom of primordial skull; B () is the image after the method process that proposes of the present invention; C () is the Yuan Tu local magnification region chosen; D () is the former figure of local magnification region; E () is the local area image after the process of adaptivelowpassfilter method; F () is the local area image after the process of wavelet-Fourierfiltering method; G () is the L that the present invention adopts ₀the result of smoothing method when Selecting parameter λ is 0.007.

Embodiment

Elaborate below in conjunction with the embodiment of the drawings and specific embodiments to the inventive method.

A kind of CBCT image annular artifact eliminating method, as shown in Figure 1, its concrete implementation step is as follows:

One, in order to change the geometrical property of annular artifact, the difficulty reduced in artifact extraction and processing procedure, image is transformed to by cartesian coordinate system the polar coordinate system that artifact shows as parallel lines on vertical direction;

For any point (x, y) in cartesian coordinate system, if the polar coordinates point of its correspondence is (ρ, θ), so their relation is:

$\left\{\begin{array}{c}{\mathrm{\ρ}}^2=x^2+y^2\\ \tan \mathrm{\θ}=\frac{y}{x}\left(x\≠0\right)\end{array}\right.---\left(1\right)$

The original-gray image F size adopted at this is 512*512's, sample conversion being carried out to image, in order to make the loss of detail of original image reach minimum in this course, smaller sampling angle and radial sampled distance must be adopted, here, we establish Δ ρ=1 and because:

Thus, the span of θ be equal to 360 ° namely 2 π time situation.Simultaneously, in order to ensure that image changes round its reconstructed center, we specify:

$\left\{\begin{array}{c}x=r\cos \mathrm{\θ}-256.5\\ y=r\sin \mathrm{\θ}-256.5\end{array}\right.---\left(3\right)$

In conjunction with formula (1), we just probably can obtain the reconstruction image P under polar coordinates.

Because P is got by sampled point conversion, perhaps, its resolution can not reach the requirement of target resolution (as 360*360), in such cases, by the most basic, interpolation point is set at the most neighbor point of known pixels point, and the most contiguous difference approach that the numerical value of interpolation point rounds up is compensated the image P under polar coordinates, just can obtain the pixel value of each point in polar coordinates accurately.

Image before and after conversion is as shown in Fig. 2 (a), (b), as seen from the figure, in polar coordinates, the annular artifact in cartesian coordinate system becomes the more simple and distinct stripe artifact of geometric properties on vertical direction, and the resolution of image there is no obvious change.

Two, L is adopted ₀the method of norm filtering obtains smoothed image P' to the smoothing process of CT image P under polar coordinate system;

According to the superiority of this kind of filtering method, it is intact and very level and smooth without artifacts P ' in certain area that we can obtain a breadths edge, only has a small amount of detailed information and be removed substantially simultaneously in this process.

Wherein L ₀norm filter achieving method is as follows:

If I is input picture, S is the Output rusults after smothing filtering, then in its gradient function, the number of nonzero value is:

C(S)＝#{p||δ _xS _p|+|δ _yS _p|≠0}(4)

Wherein p represents a pixel in S, | δ _xs _p| represent the absolute value of the gradient difference of p and x-axis direction neighbor pixel, | δ _ys _p| represent the absolute value of the gradient difference of p and y-axis direction neighbor pixel, #{} is a count operator, exports in S all satisfied | δ _xs _p|+| δ _ys _p| the number of ≠ 0 pixel.

By the L to non-zero gradient number and gradient in image ₀norm retrains, and we can in stick signal when main information, to the smoothing elimination of some unessential information.Meanwhile, similar as much as possible to original signal I in order to reach filtered signal S, we also need the reduction of gradient and structure similar between search out a balance, so there is energy function formula:

Min _s{ Σ (S _p-I _p) ²+ γ C (S) } (5) wherein Σ (S _p-I _p) ²this item constraint output image S will be similar as much as possible to input picture I, γ is a parameter for control C (S) weight, we are referred to as smoothing factor, a larger γ means in result have less edge, otherwise a less γ then means that the most detailed information in signal will be retained.

Under existing study condition, as the L of np hard problem ₀norm is unsolvable, and thus in order to obtain the optimum solution of above energy function, algorithm can be broken down into two subproblems and solve respectively, and finally obtains the more excellent solution that is similar to.

The following detailed description of once L ₀the solution procedure of norm.

Introduce two new variable h _pand v _pcorresponding δ xS respectively _pwith δ yS _p, and energy function is rewritten as:

Min _{s, h, v}{ Σ _p(S _p-I _p) ²+ γ * c (h, v)+β * ((δ xS _p-h _p) ²+ (δ yS _p-v _p) ²) (6) wherein, c (h, v)=#{p||h _p|+| v _p| ≠ 0}, β are used for controling parameters (h, v) and (δ xS as an adaptive adjustment parameter _p, δ _ys _p) similarity.In the very large situation of β, this function is equal to the energy function shown in formula (5).Because when known S, c (h, v) can separate, and when known c (h, v), S is also what can separate, as long as so we minimize c (h respectively, v) row iteration of going forward side by side with S solves, and along with β becomes large gradually, just can obtain the result of a satisfaction.

Subproblem (1): calculate S

Remove the item irrelevant with result S, then this problem equivalent is in minimizing following function:

Min _S{Σ _p(S _p-I _p) ²+β*((δxS _p-h _p) ²+(δyS _p-v _p) ²)}(7)

Utilize Fast Fourier Transform (FFT) to solve above formula and can obtain following result:

Wherein, represent Fourier transform, for inversefouriertransform, * represents conjugate function, conjugation is got after representing Fourier transform, represent the Fourier transform of impulse function, δ x and δ y represents the gradient algorithm matrix in x direction and y direction respectively, h and v be the gradient matrix δ xS of correspondence image S in x direction and the gradient matrix δ yS in y direction respectively, and adding in formula all carries out computing by matrix component one by one with multiplication and division operational symbol.Known by experiment, by this method function is solved quicker than solving function in spatial domain.

Subproblem (2): calculate (h, v)

Identical with subproblem (1), we delete the item irrelevant with (h, v), then this subproblem is equivalent to and minimizes as minor function:

${\mathrm{Min}}_{h,v}\left\{{\mathrm{\Σ}}_p\left({\left({\mathrm{\δxS}}_p-h_p\right)}^2+{\left({\mathrm{\δyS}}_p-v_p\right)}^2\right)+\frac{\mathrm{\γ}}{\mathrm{\β}}*c\left(h,v\right)\right\}---\left(9\right)$

Wherein c (h, v) represents | the number of nonzero value in h|+|v|.Due to for each independently p, h _pwith v _pseparate.So above formula be one can by the subproblem of rapid solving, we will sue for peace outside symbol extraction to formula:

${\mathrm{\Σ}}_p{\min }_{h_p,v_p}\left\{{\left(h_p-{\mathrm{\δxS}}_p\right)}^2+{\left(v_p-{\mathrm{\δyS}}_p\right)}^2+\frac{\mathrm{\γ}}{\mathrm{\β}}H\left(\left|h_p\right|+\left|v_p\right|\right)\right\}---\left(10\right)$

Wherein, H (| h _p|+| v _p|) be a two-valued function, when it returns 1 | h _p|+| v _p| ≠ 0, other then return 0.So there is sub-energy function for each pixel p:

$E_p=\left\{{\left(h_p-{\mathrm{\δxS}}_p\right)}^2+{\left(v_p-{\mathrm{\δyS}}_p\right)}^2+\frac{\mathrm{\γ}}{\mathrm{\β}}H\left(\left|h_p\right|+\left|v_p\right|\right)\right\}---\left(11\right)$

The solution of this function is:

$\left(h_p,v_p\right)=\left\{\begin{array}{cc}\left(0,0\right)& {\left({\mathrm{\δxS}}_p\right)}^2+{\left({\mathrm{\δyS}}_P\right)}^2\≤\mathrm{\γ}/\mathrm{\β}\\ \left({\mathrm{\δxS}}_p,{\mathrm{\δyS}}_P\right)& otherwise\end{array}\right.---\left(12\right)$

Thus, as long as its least energy can be calculated to each pixel so they and be exactly the globally optimal solution of formula (6).

Three, because filtering has paid the utmost attention to the removal of inessential information and stripe artifact, inevitably cause the loss of image section detailed information, in order to compensate this part of details in original image, first the smoothed image P' obtained after the original image P under polar coordinate system and filtering is subtracted each other the artifact hum pattern P obtaining having the original detailed information of part _k, then to P _kprocess, the pseudo-shadow information of image (artifact template) is further separated from detailed information, thus remaining detailed information can compensate back filtered image P ';

In order to improve separation accuracy, reduce workload, avoid as the regions of non-interest gray-scale values such as background on artifact template extraction the impact that produces, this part work is subdivided into again two stages:

1, first to the area-of-interest having obvious stripe-shaped artifact, namely as lower in gray-scale values such as soft tissues in medical image and the region with the feature such as comparatively level and smooth is within the specific limits extracted, thus reaches the object obtaining comparatively accurately artifact template.

The present invention have chosen and is considered to maximum variance between clusters (OSTU) that adaptive threshold chooses field optimum always and carries out automatic threshold value to image and choose and segmentation, obtains the area-of-interest in image.

Traditional OSTU method is a simple binary segmentation method, and for a secondary given image F, it has X different gray level, if N _ifor the number of pixels that i-th gray level has, for its total pixel number of F be so wherein the probability of i-th gray level appearance is

$\begin{array}{cc}{p}_i=\frac{N_i}{N}& p_i\≥0,{\mathrm{\Σ}}_{i=1}^X{p}_i=1---\left(13\right)\end{array}$

Namely the segmentation of single threshold is Given Graph picture is divided into object and background two large classes, sets target class as the less part of gray-scale value thus, if adaptively can choose a threshold value w, makes target class F ₀gray level [1,2 ..., w] between and background classes F ₁gray-scale value [w+1, w+2 ..., X] between, just successfully can reach the object that object and background is separated.

We provide probability and the average gray level of the appearance of each class, respectively:

${\mathrm{\ω}}_0={\mathrm{\Σ}}_{i=1}^w{p}_i=\mathrm{\ω}\left(w\right)---\left(14\right)$

${\mathrm{\ω}}_1={\mathrm{\Σ}}_{i=w+1}^X{p}_i=1-\mathrm{\ω}\left(w\right)---\left(15\right)$

And

${\mathrm{\μ}}_0=\frac{{\mathrm{\Σ}}_{i=1}^w{\mathrm{ip}}_i}{{\mathrm{\ω}}_0}=\mathrm{\μ}\left(w\right)/\mathrm{\ω}\left(w\right)---\left(16\right)$

${\mathrm{\μ}}_1=\frac{{\mathrm{\Σ}}_{i=w+1}^X{\mathrm{ip}}_i}{{\mathrm{\ω}}_1}=\frac{{\mathrm{\μ}}_{\mathrm{\τ}}-\mathrm{\μ}\left(w\right)}{1-\mathrm{\ω}\left(w\right)},{\mathrm{\μ}}_{\mathrm{\τ}}={\mathrm{\Σ}}_{i=1}^X{\mathrm{ip}}_i---\left(17\right)$

Thus, the variance of two classes is respectively:

${\mathrm{\σ}}_0^2={\mathrm{\Σ}}_{i=1}^w{\left(i-{\mathrm{\μ}}_0\right)}^2{p}_i/{\mathrm{\ω}}_0---\left(18\right)$

${\mathrm{\σ}}_1^2={\mathrm{\Σ}}_{i=w+1}^x{(i-{\mathrm{\μ}}_1)}^2{p}_i/{\mathrm{\ω}}_1---\left(19\right)$

Regulation maximum between-cluster variance is:

σ ²=ω ₀(μ ₀-μ _τ) ²+ ω ₁(μ ₁-μ _τ) ²=ω ₀ω ₁(μ ₁-μ ₀) ²(20) then for each w, have:

${\mathrm{\σ}}^2\left(w\right)=\frac{{\[{\mathrm{\μ}}_{\mathrm{\τ}}\mathrm{\ω}\left(w\right)-\mathrm{\μ}\left(w\right)\]}^2}{\mathrm{\ω}\left(w\right)\[1-\mathrm{\ω}\left(w\right)\]}---\left(21\right)$

Wherein optimum w ^*for:

${\mathrm{\σ}}^2\left(w^{*}\right)=\underset{1\≤w\≤X}{\max }{\mathrm{\σ}}^2\left(w\right)---\left(22\right)$

In real medical image, except background, some specific region is as human skeleton, its gray-scale value will have very large lifting compared to regions such as soft tissues, in these regions, we often do not observe the existence of annular artifact, and the violent conversion of edge's gray-scale value can produce certain impact to the generation of artifact template on the contrary, regions of non-interest also will be thought by us in these regions thus, in order to eliminate this impact, we will carry out multi-threshold segmentation to image further, further piecemeal process is carried out to two large classes in former methodical basis thus reaches the object of image being carried out to multi-threshold segmentation.For image P, if it will be divided into m region, so we will just need m-1 threshold value to divide it.According to character and the many experiments of tissue conical beam CT image, we select m=3.In concrete operating process, for already present class, we find the maximum local of variance within clusters in the overall situation and carry out the group of OSTU Threshold segmentation as the next one and newly-generated class recalculated to its variance within clusters, by that analogy until obtain m-1 threshold value needed for us.The computing formula of the variance within clusters of local is:

${\mathrm{\σ}}^2={\mathrm{\ω}}_0{\mathrm{\σ}}_0^2+{\mathrm{\ω}}_1{\mathrm{\σ}}_1^2---\left(23\right)$

Further, the optimal threshold that this group obtained contains m-1 numeral still meets the maximum principle of inter-class variance, that is:

σ ²＝ω ₀(μ ₀-μ _τ) ²+ω ₁(μ ₁-μ _τ) ²+…+ω _m-1(μ _m-1-μ _τ) ²(24)

By this kind of method, we just successfully by needing the tissue region of carrying out the extraction of stripe-shaped artifact to separate from the background of image, both can improve the accuracy of artifact template generation, having turn eliminated the impact that background area pixels value may be brought.Because OSTU algorithm is inherently a kind of exhaustive searching algorithm; so in the writing process of program; before the OSTU segmentation carrying out local; usually Two-peak method can be utilized to determine the scope of a gray-scale value as hunting zone during threshold search to reduce the time of search, and we claim this scope for lax surplus.

It should be noted that, find after many sub-pictures are tested, multi-threshold segmentation is done on the original image P in polar coordinates can not eliminate with process the impact that regions of non-interest causes template extraction completely (this impact is shown to be due to L through experimental analysis only separately ₀norm filtering some limitation on edge and background process caused).Thus, in order to obtain the position of area-of-interest in conical beam CT image accurately, eliminating the impact in the region such as background and bone, generating artifact template more really, we are to the P of former figure P and artifact hum pattern _kdo multi-threshold segmentation respectively and in area-of-interest, get it and occur simultaneously, show in the area-of-interest extracted in this mode through test of many times, the latter will obtain the most accurate artifact template.

2, based on area-of-interest, at the P of artifact hum pattern _kthis region in utilize median method to carry out extracting and remove to artifact template.

From annular artifact feature, for the different pixels point on annular artifact caused by same detector, its gray-scale value is substantially identical.Thus, we intermediate value is got for pixel of each perpendicular row in area-of-interest and generate thus one with the artifact template of area-of-interest with size, thereafter, with the P of artifact hum pattern _kdeduct artifact template and just can obtain filtered affected detailed information in image, thus compensated back in filtered image P'.

In practical operation, in order to get rid of the impact (so-called bad pixel is the pixel be adjacent and compares the too violent pixel of gray-value variation) that some bad pixels cause, in the process of getting intermediate value, get rid of higher with the lower point of some gray-scale values and often obtain better effect (being generally front and back 2%).

Four, the image P ' that compensate for detailed information is switched back to cartesian coordinate system, obtain the conical beam CT image processing rear ringless-type artifact.

Type in polar domain is double and image P ' (ρ, θ) after the process of grey scale pixel value all between 0 to 1, if its size is M*N (in the image be applied to herein, M=N=360), according to formula:

$\left\{\begin{array}{c}x=\mathrm{\ρ}\cos \mathrm{\θ}\\ y=\mathrm{\ρ}\sin \mathrm{\θ}\end{array}\right.---\left(25\right)$

Convert it back to cartesian coordinate system, it should be noted that in cartesian coordinate system, the size of image, with originally identical, is 512*512.In this step, we carry out Interpolation compensation process by utilizing the method for bilinear interpolation to image thus, finally obtain ringless-type artifact under cartesian coordinate system and image resolution ratio rebuild without the high-quality obviously declined after medical science conical beam CT image.

Embodiment:

Conical beam CT image totally four width that experiment adopts, be the real human body tissue of utilization captured by operation CT machine, be respectively the basis cranii CT image of two width Cranial Computed Tomography images, a width neck CT image and a width downwards angle of visibility, these images are by part compositions such as a large amount of human body soft tissues, there will be comparatively significantly annular artifact, is the keypoint part of current industry research.Conveniently process, four width images are the gray level image of 512*512 size, and at L ₀in the process of norm filtering, β ₀with β _maxunified employing 2 λ and 1E5, iterative rate gets 2, and filtering parameter λ then needs to adjust within the specific limits according to filter effect.From Fig. 3,4, identical position can select identical parameter carry out processing and obtain good effect with the CT image of angle; And for the image of different angles different parts, filtering parameter then needs to be determined by experiment.Empirical tests, this scope is usually between 0.005 to 0.03.

Experimental result as depicted in figs. 3-6, the algorithm we adopted is with the comparatively classical annular artifact minimizing technology (wavelet-Fourierfiltering) based on wavelet decomposition and Fourier transform and compare based on the method (adaptivelowpassfilter) of unidirectional adaptive smooth process, with regard to integral image, of the present invention by L ₀the annular artifact minimizing technology that norm filtering combines with artifact template extraction to retain in conical beam CT image while details and marginal information to the full extent, effectively eliminates the annular artifact that it affects picture quality and authenticity and there is no obviously residual.Meanwhile, except the overall gray value of image there occurs except slight conversion, there is not obvious decline in the resolution of image itself yet.It should be noted that the slight problem be not both caused by the method for interpolation in coordinate transformation process but not existing for method itself of some fringe regions.

Above technical scheme of the present invention is described; but these explanations can not be understood to limit scope of the present invention; protection scope of the present invention is limited by the claims of enclosing, and any change on the claims in the present invention basis is all protection scope of the present invention.

Claims (5)

1. a CBCT image annular artifact eliminating method, is characterized in that: comprise the following steps:

Step one, transforms to artifact by original image I by cartesian coordinate system and shows as the image P that the polar coordinate system of parallel lines on vertical direction obtains under polar coordinate system;

Step 2, adopts L ₀the method of norm filtering obtains smoothed image P' to the smoothing process of image P under polar coordinate system;

Step 3, by subtracting each other the artifact hum pattern P obtaining having the original detailed information of part by the smoothed image P' obtained after the original image P under polar coordinate system and filtering _k, then to P _kextraction and the removal of carrying out artifact template obtain P _k', finally will eliminate artifact template and only comprise the P of detailed information _k' compensate back filtered image P ';

Step 4, switches back to cartesian coordinate system by the image P ' that compensate for detailed information, just can obtain the cone beam images of a width ringless-type artifact.

2. a kind of CBCT image annular artifact eliminating method according to claim 1, is characterized in that: after obtaining image P described in step one, adopt the most contiguous difference approach to compensate to reach resolution requirement to it.

3. a kind of CBCT image annular artifact eliminating method according to claim 1 and 2, is characterized in that: to P described in step 3 _kcarry out process obtain artifact template and it removal is completed by following process:

(1) OSTU multi-threshold segmentation method is adopted to extract the area-of-interest having obvious stripe-shaped artifact;

(2) at artifact hum pattern P _karea-of-interest in utilize median method to carry out extracting and remove to artifact template.

4. a kind of CBCT image annular artifact eliminating method according to claim 3, is characterized in that: the described area-of-interest to having obvious stripe-shaped artifact of step (1) is extracted and completed by following process: to former figure P and artifact hum pattern P _kdo multi-threshold segmentation respectively and in area-of-interest, get its common factor and obtain final area-of-interest.

5. a kind of CBCT image annular artifact eliminating method according to claim 3, is characterized in that: the described median method that utilizes of step (2) carries out extracting to artifact template and removes further by getting the front-seat point higher with lower except some gray-scale values of intermediate value to obtain better effect.