CN102142133A - Mammary X-ray image enhancement method based on non-subsampled Directionlet transform and compressive sensing - Google Patents

Abstract

The invention discloses a mammary X-ray image enhancement method based on non-subsampled Directionlet transform and compressive sensing, and the method is mainly used for solving the defect of insufficient enhancement effect on the low-contrast medical image in the existing method. The image enhancement method is characterized by comprising the following steps: introducing non-subsampled Directionlet transform and compressive sensing into image enhancement, namely, firstly performing non-subsampled Directionlet transform on the image, and centralizing energy with a high-frequency coefficient by utilizing a compressive sensing technology; then enhancing the concentrated high-frequency coefficient by utilizing a linear enhancement algorithm; and finally reconstructing the enhanced frequency domain representation coefficient through non-subsampled Directionlet transform so as to obtain the enhanced mammary image. By utilizing the mammary X-ray image enhancement method, influence of pathological change region background can be better inhibited, features of a pathological change region in the low-contrast image are obviously enhanced, and the information quantity and readability of the image are improved, thus the method can be used in aided medical radiodiagnosis.

Description

Mammary X-ray image enchancing method based on non-lower sampling Directionlet conversion and compressed sensing

Technical field:

The invention belongs to image processing field, relate to the image enchancing method of non-lower sampling Directionlet conversion and compressed sensing, this method can strengthen the mammary X-ray image of low contrast effectively, and the doctor of auxiliary radiation section carries out medical diagnosis.

Background technology:

Breast cancer is comparatively one of common malignancy disease of women, and women's life and health in serious threat.Development along with the modern medicine imaging technique, various new medical imaging technologies have been widely used in each links such as medical diagnosis, the preceding plan of art, treatment, monitoring after operation, these imaging techniques can obtain patient's various data comprehensively and exactly, for diagnosis, treatment, operation and postoperative evaluation provide information more accurately, wherein the soft X line of molybdenum target is higher with its spatial resolution, comparatively responsive to lump and calcification, and advantage such as equipment needed thereby is simple, cheap, become present early-stage breast cancer diagnosis means the most reliable and commonly used.Yet misdiagnosis rate and rate of missed diagnosis when using the soft X-ray film of breast molybdenum target to diagnose are still higher, and this mainly is because inferior picture quality, the optimum performance of malignant change and observer's visual fatigue or carelessness.Wherein inferior picture quality is embodied in the following aspects: the area-of-interest contrast is relatively poor, suspicious lesions is regional and its surrounding tissue between intensity difference very faint, the lesion region shape is changeable and not of uniform size and obscurity boundary etc.Along with computing machine and development of technologies thereof, utilize computer technology that the mammary X-ray image is strengthened and can effectively address this problem, the help doctor understands image better and judges.

In order to highlight the feature of lesion region, improve the visual effect of mammary X-ray image, the most frequently used method is carried out enhancement process to image exactly.Traditional image enhancement processing technology has obtained to a certain extent and has strengthened effect preferably, but for the mammary X-ray image that has than low contrast, their enhancing effect can not be satisfactory.So far, be directed to the mammary X-ray image, proposed multiple Enhancement Method.

1. anti-sharpening masking method

Anti-sharpening masking method is one of Enhancement Method commonly used in the Flame Image Process, and it is meant on the basis of original image, adds a certain proportion of image radio-frequency component, in the hope of reaching the effect that edge and detailed information strengthen.The advantage of anti-sharpening masking method is to give prominence to edge of image and detailed information preferably, thereby reaches the enhancing effect to image.But this anti-sharpening masking method is owing to the contrast information that does not take into full account image, thereby under the lower situation of contrast, it is unsatisfactory that it strengthens effect.

2. adaptive histogram equalization method

The adaptive histogram equalization method is the effective image enchancing methods of a kind of classics, it adopts sliding window technique, the window area that comprises processed point is carried out histogram equalization, be about to the region histogram distribution and change into even histogram distribution, and give pixel pending in window assignment again according to the mapping relations of histogram and gray scale on this basis.The advantage of adaptive histogram equalization method is the dynamic range that can be good at regulating image, strengthens the details of image simultaneously, but this method has also been amplified noise when improving contrast, so its enhancing effect is still waiting to improve.

3. wavelet transformation adaptive gain disposal route

Wavelet transformation adaptive gain disposal route is a kind of image enchancing method comparatively commonly used, it at first carries out wavelet transformation to image, then according to the adaptive enhancement process of carrying out of the distribution situation of image transformation coefficient, again the conversion coefficient after handling is carried out inverse transformation, the image after being enhanced at last.The advantage of wavelet transformation adaptive gain disposal route is to improve the contrast of image preferably, and noise is had certain robustness, but when the gray difference of the lesion region of image and background area was very little, its enhancing effect still had much room for improvement.

Though three kinds of above-mentioned methods have obtained to a certain extent and have strengthened effect preferably, but because the mammary X-ray image has that contrast is low, noise is than characteristics such as the gray difference of horn of plenty and lesion region and background area are little, so said method is not very good to the enhancing effect of mammary X-ray image.

Summary of the invention

The objective of the invention is to overcome the deficiency of above-mentioned prior art, a kind of mammary X-ray image enchancing method based on non-lower sampling Directionlet conversion and compressed sensing has been proposed, background with effective inhibition image, highlight the lesion region feature, improve the contrast of image, make the enhancing effect of mammary X-ray image more obvious.

Realize that the object of the invention technical thought is: by removing the down-sampling operation in the Directionlet conversion, realized the Directionlet conversion of non-lower sampling, and it is combined with compression sensing method, strengthened clear effect to the mammary X-ray image.Its concrete scheme comprises the steps:

(1) with input picture I _InTransform to the Directionlet transform domain, promptly utilize non-lower sampling Directionlet transfer pair input picture I _InCarry out sub-band division, obtain frequency domain representation coefficient D, this frequency domain representation coefficient D comprises high fdrequency component

And low frequency component

(2) generate white Gaussian noise observing matrix Φ at random.

(3) use the white Gaussian noise observing matrix Φ of generation to the high fdrequency component among the frequency domain representation coefficient D

Observe, obtain observed reading X.

(4) adopt the OMP algorithm, X recovers to observed reading, the high fdrequency component after being restored

(5) to the high fdrequency component after recovering Carry out linear enhancement process, the high fdrequency component after being enhanced

(6) to the low frequency component among the frequency domain representation coefficient D And the high fdrequency component after strengthening

Carry out non-lower sampling Directionlet inverse transformation, the image I after finally being enhanced _Out

The present invention has the following advantages:

(1) the present invention has been owing to adopted the Directionlet transform method, thereby can catch the direction detailed information of image effectively, highlighted the minutia of lesion region in the mammary X-ray image significantly, improves the sharpness of image.

(2) the present invention concentrates image energy effectively, thereby can carry out enhancement process to high-frequency information in the concentrated area owing to adopt high fdrequency component after the method for compressed sensing is restored, improves the contrast of image significantly.

(3) the present invention is simple in structure, and calculation cost is low, can carry out enhancement process to image simply and effectively.

Description of drawings:

Fig. 1 is a process flow diagram of the present invention.

Fig. 2 is the schematic diagram that the present invention uses sampling matrix that image is sampled.

Fig. 3 is that the present invention uses the figure as a result after sampling matrix is sampled to image.

Fig. 4 is that the present invention uses sampling matrix frequency coefficient to be carried out the sub-process figure of matrix interpolation and stack.

Fig. 5 uses the present invention and existing method to mammary X-ray treatment of picture effect contrast figure.

Specific embodiments:

With reference to Fig. 1, the present invention includes non-lower sampling Directionlet conversion, compressed sensing, enhancement process and non-lower sampling Directionlet inverse transformation.Concrete steps are as follows:

Step 1: generate sampling matrix.

1.1) be starting point with the true origin, produce the integer vectors d of two linear independences at random ₁And d ₂, respectively as image changing direction and formation direction, wherein d ₁=[a ₁, b ₁], d ₂=[a ₂, b ₂], a ₁Be integer vectors d ₁The terminal point horizontal ordinate, b ₁Be integer vectors d ₁The terminal point ordinate, a ₂Be integer vectors d ₂The terminal point horizontal ordinate, b ₂Be integer vectors d ₂The terminal point ordinate;

1.2) with integer vectors d ₁And d ₂Constitute sampling matrix M _Λ:

$M_{\mathrm{\Λ}}=\left(\begin{array}{c}{d}_1\\ d_2\end{array}\right)=\left(\begin{array}{cc}{a}_1& b_1\\ a_2& b_2\end{array}\right),a_1,a_2,b_1,b_2\⋐Z$

In the formula, Z is the integer complete or collected works.

Step 2: according to sampling matrix M _ΛWith input picture I _InBe divided into | det (M _Λ) | individual mutual incoherent subgraph sequence F, each subgraph among the F is F _k, k=0 ..., | det (M _Λ) |-1,

Wherein, | det (M _Λ) | be sampling matrix M _ΛDeterminant, F _kBe arbitrary subgraph, be expressed as: F _k(n)=I _In(M _Λ(n-S _k))

In the formula, n=(n ₁n ₂) be the position coordinates of pixel in image, F _k(n) the position coordinate is (n in k subgraph of expression ₁n ₂) pixel, S _kRepresent k subgraph F _kDisplacement vector, I _In(M _Λ(n-S _k)) be k subgraph, it is to utilize displacement vector S _kWith input picture I _InPixel be shifted, re-use sampling matrix M _ΛIt is carried out obtaining after the matrix sampling computing, and sampling principle is chosen integer vectors d respectively as shown in Figure 2 among Fig. 2 ₁=[1,1] and d ₁=[1,1] is as changing direction and the formation direction, and the sampling matrix of generation is:

Utilize displacement vector S _kTo input picture I _InPosition coordinates be shifted after, re-use the sampling matrix M of generation _ΛWith the displacement after the picture position coordinate multiply each other, the subgraph sequence after obtaining sampling as shown in Figure 3, sampling matrix M _ΛInput picture is divided into two subgraphs, and wherein stain is represented subgraph F ₁, white point is represented subgraph F ₂

Step 3: antithetical phrase graphic sequence F carries out three layers of redundant wavelet decomposition, obtains the high-frequency sub-band coefficient sequence of subgraph sequence F under each yardstick: W ^j=(H ^j, V ^j, D ^j), j=1,2,3 and low frequency sub-band coefficient sequence A ³, high-frequency sub-band coefficient sequence W ^jIn each sub-band coefficients be:

K=0 ..., | det (M _Λ) |-1, low frequency sub-band coefficient sequence A ³In each sub-band coefficients be:

K=0 ..., | det) M _Λ) |-1;

Wherein,

Be meant expression subgraph F after carrying out the redundant wavelet decomposition of j layer _kThe frequency domain sub-band coefficients of middle horizontal direction detailed information,

Be meant expression subgraph F after carrying out the redundant wavelet decomposition of j layer _kThe frequency domain sub-band coefficients of middle vertical direction detailed information,

Be meant expression subgraph F after carrying out the redundant wavelet decomposition of j layer _kThe frequency domain sub-band coefficients of middle diagonal detailed information,

Be meant expression subgraph F after carrying out the 3rd layer of redundant wavelet decomposition _kThe frequency domain sub-band coefficients of medium and low frequency profile information, frequency domain sub-band coefficients under the j yardstick Be expressed as successively:

$H_k^j=A_k^{j-1}*h_1^{j-1}*h_0^{j-1}$

$V_k^j=A_k^{j-1}*h_0^{j-1}*h_1^{j-1}$

$D_k^j=A_k^{j-1}*h_1^{j-1}*h_1^{j-1}$

$A_k^j=A_k^{j-1}*h_0^{j-1}*h_0^{j-1},$ j＝1，2，3

In the formula,

Be the low pass resolution filter coefficient of " 97 " small echo, Be the high pass resolution filter coefficient of " 97 " small echo,

With Be respectively right

With

Carry out the filter coefficient that j interpolation obtains,

Represent subgraph to be transformed.

Step 4: according to sampling matrix M _ΛRespectively to high-frequency sub-band coefficient sequence W ^j=(H ^j, V ^j, D ^j) and low frequency sub-band coefficient sequence A ³By matrix interpolation and stack, obtain input picture I _InCarry out the frequency domain representation coefficient D=(DH after the non-lower sampling Directionlet conversion ^j, DV ^j, DD ^j, DA ³).

Referring to Fig. 4, being implemented as follows of this step:

4.1) the matrix interpolation computing

Respectively to high-frequency sub-band coefficient sequence W ^j=(H ^j, V ^j, D ^j) and low frequency sub-band coefficient sequence A ³Carrying out sampling matrix is M _ΛThe matrix interpolation computing, obtain the high-frequency sub-band coefficient sequence after the interpolation: W ^{J '}=(H ^{J '}, V ^{J '}, D ^{J '}) and low frequency sub-band coefficient sequence A ^{3 '}, be expressed as:

j＝1，2，3

In the formula, M _ΛBe sampling matrix, n=(n ₁n ₂) position coordinates of remarked pixel point in image,

Expression is carried out the new position coordinates that obtains after the matrix interpolation computing to the position coordinates of pixel, and Λ represents the subscript collection of frequency domain sub-band coefficients;

4.2) stack of sub-band coefficients sequence

Respectively to the high-frequency sub-band coefficient sequence W behind the interpolation arithmetic ^{J '}=(H ^{J '}, V ^{J '}, D ^{J '}) and low frequency sub-band coefficient sequence A ^{3 '}Superpose, obtain input picture I _InCarry out the frequency coefficient D=(DH after the non-lower sampling Directionlet conversion ^j, DV ^j, DD ^j, DA ³),

Wherein: DH ^jBe meant the change direction high fdrequency component of detailed information of presentation video, DV ^jBe meant the high fdrequency component of presentation video formation direction detailed information, DD ^jBe meant that presentation video is changed direction and the high fdrequency component of formation direction detailed information, DA ³Be meant the low frequency component of presentation video low frequency profile information, they are expressed as respectively:

${\mathrm{DH}}^j\left(n\right)=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\Λ}}\right)\right|-1}{\mathrm{\Σ}}}{H}_k^{j^{\′}}(n+S_k),$

${\mathrm{DV}}^j\left(n\right)=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\Λ}}\right)\right|-1}{\mathrm{\Σ}}}{V}_k^{j^{\′}}(n+S_k),$

${\mathrm{DD}}^j\left(n\right)=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\Λ}}\right)\right|-1}{\mathrm{\Σ}}}{D}_k^{j^{\′}}(n+S_k),$

${\mathrm{DA}}^3\left(n\right)=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\Λ}}\right)\right|-1}{\mathrm{\Σ}}}{A}_k^{3^{\′}}(n+S_k),$ j＝1，2，3

In the formula, n=(n ₁n ₂) position coordinates of remarked pixel point in image, S _kRepresent k displacement vector, (n+S _k) expression new position coordinates that position coordinates is shifted and obtains, DH ^j, DV ^j, DD ^jForm the high fdrequency component among the frequency domain representation coefficient D DA ³Form the low frequency component among the frequency domain representation coefficient D

Step 5: generate a gaussian random matrix Φ at random as observing matrix, the size of Φ is M * N, wherein M is an observation frequency, N is the line number that needs the frequency domain representation coefficient of observation, the present invention is on the basis of carrying out repeatedly experiment test and checking, the span that draws observation frequency M is between 250～500, and the observation frequency of this example is taken as 400.

Step 6: use observing matrix Φ that the high fdrequency component among the frequency domain representation coefficient D is observed, promptly with observing matrix Φ that generates at random and the high fdrequency component among the frequency domain representation coefficient D

Multiply each other, obtain observed reading

Wherein:

Be to high fdrequency component

In the change direction coefficient DH of detailed information of presentation video ^jThe result who observes is expressed as:

Be to high fdrequency component

The coefficient DV of middle presentation video formation direction detailed information ^jThe result who observes is expressed as:

Be to high fdrequency component

In presentation video change direction and the coefficient DD of formation direction detailed information ^jThe result who observes is expressed as:

Step 7: the observing matrix Φ that utilizes step 2 to generate, use orthogonal matching pursuit OMP method that observed reading X is recovered, the high frequency coefficient after the recovery is

7.1) iteration count t is initialized as 1, i.e. t=1;

7.2) the calculating observation matrix

In the column vector of each row

With observed reading

In i part coefficient

P row column vector v _pBetween inner product, obtain inner product sequence I, each inner product is I among the I _jBe expressed as:

p＝1，2，...，n，j＝1，2，...，N

In the formula, n is a coefficient

Columns, N is the columns of observing matrix Φ;

7.3) calculate the inner product of numerical value maximum among the inner product sequence I:

7.4) calculating maximum inner product I _MaxPosition number in inner product sequence I:

7.5) the maximum inner product I of usefulness _MaxPosition number λ _tConstitute sequence number collection: Λ _t=(Λ _T-1, λ _t), in the formula, sequence number collection Λ ₀Be empty set:

The expression empty set;

7.6) with λ among the observing matrix Φ _tThe column vector of row

Constitute augmented matrix:

In the formula, augmented matrix Φ ₀Be empty set:

7.7) with column vector

From observing matrix Φ, remove, be about to λ among the observing matrix Φ _tThe column vector of row

Be changed to empty set:

7.8) utilize augmented matrix Φ _tTo column vector v _pEstimate, obtain new estimated value x _t:

$x_t=\frac{<{\mathrm{\&Phi;}}_t,v_p>}{<{\mathrm{\&Phi;}}_t,{\mathrm{\&Phi;}}_t>}$

In the formula,＜Φ _t, v _pBe meant augmented matrix Φ _tWith column vector v _pInner product,＜Φ _t, Φ _tBe meant augmented matrix Φ _tWith the inner product of itself;

7.9) according to new estimated value x _tWith position number λ _tCalculate the column vector after reconstruct recovers:

7.10) utilize new estimated value x _t, augmented matrix Φ _tWith column vector v _pCalculate residual values r _t: r _t=v _p-Φ _tx _t

7.11) increase iteration count t:t=t+1, use residual values r _tReplacement step 7.2) and 7.8) in column vector v _p, repeating step 7.2)～7.11) equal 1/4th of observation frequency M up to iteration count t, promptly

Till, the column vector that obtain this moment

P=1,2 ..., n has constituted through the high frequency coefficient after the reconstruct recovery

Be expressed as

N is a coefficient

Columns, the high frequency coefficient after reconstruct recovers

Constituted the frequency domain high fdrequency component:

This frequency domain high fdrequency component is the main high fdrequency component feature of representative image more concentratedly, and then improves the enhancing effect of subsequent treatment to image.

Step 8: to the high frequency coefficient after the reconstruct recovery

Multiply by amplification factor μ, the high frequency coefficient after being enhanced is: Be expressed as

The present invention is on the basis of carrying out repeatedly experiment test and checking, and the span that draws reinforcing coefficient μ is between 3～6, and the observation frequency of this example is taken as 4.

Step 9: use sampling matrix M _ΛRespectively to the high fdrequency component after strengthening

With the low frequency component among the frequency domain representation coefficient D

Carry out the matrix sampling, obtain the frequency coefficient sequence

Wherein:

Be to use sampling matrix M _ΛTo high fdrequency component

Carry out the result of matrix sampling, Be to use sampling matrix M _ΛTo low frequency component DA ³Carry out the result of matrix sampling, be expressed as

${\hat{Z}}_{\mathrm{ik}}^j\left(n\right)={\hat{Y}}_i^j\left(M_{\mathrm{\&Lambda;}}(n-S_k)\right),$ i＝1，2，3

$Z_k^3\left(n\right)={\mathrm{DA}}^3\left(M_{\mathrm{\&Lambda;}}(n-S_k)\right),$ k＝0，...，|det(M _Λ)|-1

In the formula, n=(n ₁n ₂) be the position coordinates of pixel in image, S _kRepresent k displacement vector, M _Λ(n-S _k) represent by displacement vector S _kPosition coordinates is shifted, and uses sampling matrix the position coordinates after being shifted to be carried out the new position coordinates that obtains after the matrix sampling.

Step 10: with the frequency coefficient sequence

Carry out three layers of redundant wavelet reconstruction respectively, obtain the coefficient sequence Z behind the redundant wavelet reconstruction of j layer ^j, coefficient sequence Z ^jIn each subsystem number

Be expressed as:

$Z_k^j=Z_k^{j+1}*g_0^j*g_0^j+{\hat{Z}}_{3k}^{j+1}*g_1^j*g_1^j$

$+{\hat{Z}}_{2k}^{j+1}*g_1^j*g_0^j+{\hat{Z}}_{1k}^{j+1}*g_0^j*g_0^j,$ j＝2，1，0

In the formula,

With Reconstruct low-pass filter coefficients and the reconstruct Hi-pass filter coefficient of representing " 97 " small echo respectively,

With Be respectively right

With

Carry out the filter coefficient that j interpolation obtains; Obtain the image sequence after the conversion: Z this moment ⁰=Z ^j(j=0).

Step 11: according to sampling matrix M _ΛRespectively to the image sequence Z after the conversion ⁰In each subgraph

By matrix interpolation and stack, obtain through the enhancing image I after non-lower sampling Directionlet conversion and the compressed sensing enhancing _Out

11.1) matrix interpolation

To image sequence Z ⁰In each subgraph

Carrying out sampling matrix is M _ΛThe matrix interpolation computing, obtain the image sequence Z after the interpolation, each subgraph Z among the image sequence Z _kBe expressed as:

In the formula, M _ΛBe sampling matrix, n=(n ₁n ₂) position coordinates of remarked pixel point in image,

The new position coordinates that obtains after the matrix interpolation is carried out in expression to position coordinates;

11.2) the image sequence stack

Image sequence Z after the matrix interpolation computing is superposeed, obtain through the enhancing image I after non-lower sampling Directionlet conversion and the compressed sensing enhancing _Out, be expressed as:

$I_{\mathrm{out}}=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\&Lambda;}}\right)\right|-1}{\mathrm{\&Sigma;}}}{Z}_k(n+S_k)$

In the formula, n=(n ₁n ₂) position coordinates of remarked pixel point in image, S _kBe meant k displacement vector, Z _k(n+S _k) represent image sequence Z _kIn k subgraph carry out the new subgraph that obtains after the coordinate displacement.

Advantage of the present invention can further specify by following emulation experiment:

1. simulated conditions

Test pattern used in the present invention comes from the mammary X-ray image in the MIAS database.

2. emulation content

2.1) the present invention uses the sharpening mask method respectively, the adaptive histogram equalization method, these four kinds of methods of wavelet transformation adaptive gain facture and the inventive method have carried out strengthening test to four groups of mammary X-ray images, obtain four groups and strengthen the back image, as shown in Figure 5, wherein Fig. 5 (a) from top to bottom is four groups of original images successively, Fig. 5 (b) is to use existing sharpening mask method to the enhancing of Fig. 5 (a) figure as a result from top to bottom successively, Fig. 5 (c) is to use existing adaptive histogram equalization method to the enhancing of Fig. 5 (a) figure as a result from top to bottom successively, Fig. 5 (d) is to use existing wavelet transformation adaptive gain disposal route to the enhancing of Fig. 5 (a) figure as a result from top to bottom successively, and Fig. 5 (e) is to use the present invention to the enhancing of Fig. 5 (a) figure as a result from top to bottom successively.As seen Fig. 5 (e) is compared with Fig. 5 (b), Fig. 5 (c), Fig. 5 (d) respectively, the present invention can be effectively strengthens the lesion region of galactophore image, suppressed to belong in the image normal structure of background preferably for the influence that strengthens lesion region, make the complex background zone become smoothly, thereby further increased information content of image and readability.

2.2) the present invention is with detection target and background contrast ratio TB based on variance _cValue to using four groups of mammary X-ray images after these four kinds of methods of sharpening mask method, adaptive histogram equalization method, wavelet transformation adaptive gain facture and the inventive method strengthen to test, obtains four groups of contrast ratio TB respectively as judging basis _cValue, as shown in table 1, wherein first row are four groups of original sequence, TB after secondary series is to use the sharpening mask method to four picture group image intensifyings _cValue, TB after the 3rd row are to use the adaptive histogram equalization method to four picture group image intensifyings _cValue, TB after the 4th row are to use wavelet transformation adaptive gain facture to four picture group image intensifyings _cValue, TB after the 5th row are to use the inventive method to four picture group image intensifyings _cValue, wherein, contrast ratio TB _cValue representation is:

TB _c＝δ _μ/σ

In the formula, δ _μMeasured the difference between the average gray ratio that detects target and background in original image and the enhancing back image, δ _μBe expressed as: Wherein

With

Be meant the average that detects target T and background B,

With Be meant the average that strengthens the back image.σ has measured to strengthen and has detected target in the image with respect to the reduction degree that detects target gray level divergence in the original image,

In the formula

With

Be respectively original image and strengthen the variance that detects target in the image of back, therefore detect target and background contrast ratio TB _cValue is big more to show that then the enhancing effect of image is good more.

Table 1 strengthens evaluation of result TB _c

By table 1 as seen, the TB of the present invention after four groups of images are tested _cValue is obviously greater than other three kinds of existing methods, and therefore enhancing effect of the present invention is better than these three kinds of existing methods of sharpening mask method, adaptive histogram equalization method and wavelet transformation adaptive gain disposal route on objective metric.

To sum up, the present invention has improved the contrast that detects target and background, has strengthened the detailed information of galactophore image effectively.

Claims (6)

1. the mammary X-ray image enchancing method based on non-lower sampling Directionlet conversion and compressed sensing comprises the steps:

(1) with input picture I _InTransform to the Directionlet transform domain, promptly utilize non-lower sampling Directionlet transfer pair input picture I _InCarry out sub-band division, obtain frequency domain representation coefficient D, this frequency domain representation coefficient D comprises high fdrequency component

And low frequency component

(2) generate white Gaussian noise observing matrix Φ at random.

(3) use the white Gaussian noise observing matrix Φ of generation to the high fdrequency component among the frequency domain representation coefficient D

Observe, obtain observed reading X.

(4) adopt the OMP algorithm, X recovers to observed reading, the high fdrequency component after being restored

(5) to the high fdrequency component after recovering

Carry out linear enhancement process, the high fdrequency component after being enhanced

(6) to the low frequency component among the frequency domain representation coefficient D

And the high fdrequency component after strengthening

Carry out non-lower sampling Directionlet inverse transformation, the image I after finally being enhanced _Out

2. according to claims 1 described mammary X-ray image enchancing method, the described non-lower sampling Directionlet transfer pair input picture I that utilizes of step (1) wherein _InCarry out sub-band division, carry out as follows:

(1a) produce the integer vectors d of two linear independences at random ₁And d ₂, respectively as image changing direction and formation direction, wherein d ₁=[a ₁, b ₁], d ₂=[a ₂, b ₂], a ₁Be integer vectors d ₁The terminal point horizontal ordinate, b ₁Be integer vectors d ₁The terminal point ordinate, a ₂Be integer vectors d ₂The terminal point horizontal ordinate, b ₂Be integer vectors d ₂The terminal point ordinate;

(1b) with integer vectors d ₁And d ₂Constitute sampling matrix M _Λ:

$M_{\mathrm{\&Lambda;}}=\left(\begin{array}{c}{d}_1\\ d_2\end{array}\right)=\left(\begin{array}{cc}{a}_1& b_1\\ a_2& b_2\end{array}\right),a_1,a_2,b_1,b_2\&Subset;Z$

In the formula, Z is the integer complete or collected works;

(1c) according to sampling matrix M _ΛWith input picture I _InBe divided into | det (M _Λ) | individual mutual incoherent subgraph sequence F, each subgraph among the F is F _k, k=0 ..., | det (M _Λ) |-1,

Wherein, | det (M _Λ) | be sampling matrix M _ΛDeterminant, F _kBe arbitrary subgraph, be expressed as:

F _k(n)＝I _in(M _Λ(n-S _k))

In the formula, n=(n ₁n ₂) be the position coordinates of pixel in image, F _k(n) the position coordinate is (n in k subgraph of expression ₁n ₂) pixel, S _kRepresent k subgraph F _kDisplacement vector, I _In(M _Λ(n-S _k)) be k subgraph, it is to utilize displacement vector S _kWith input picture I _InPixel be shifted, re-use sampling matrix M _ΛIt is carried out obtaining after the matrix sampling computing;

(1d) antithetical phrase graphic sequence F carries out three layers of redundant wavelet decomposition, obtains the high-frequency sub-band coefficient sequence of subgraph sequence F under each yardstick: W ^j=(H ^j, V ^j, D ^j), j=1,2,3 and low frequency sub-band coefficient sequence A ³, high-frequency sub-band coefficient sequence W ^jIn each sub-band coefficients be:

K=0 ..., | det (M _Λ) |-1, low frequency sub-band coefficient sequence A ³In each sub-band coefficients be:

K=0 ..., | det (M _Λ) |-1;

Wherein,

Be meant expression subgraph F after carrying out the redundant wavelet decomposition of j layer _kThe frequency domain sub-band coefficients of middle horizontal direction detailed information,

Be meant expression subgraph F after carrying out the redundant wavelet decomposition of j layer _kThe frequency domain sub-band coefficients of middle vertical direction detailed information,

Be meant expression subgraph F after carrying out the redundant wavelet decomposition of j layer _kThe frequency domain sub-band coefficients of middle diagonal detailed information,

Be meant expression subgraph F after carrying out the 3rd layer of redundant wavelet decomposition _kThe frequency domain sub-band coefficients of medium and low frequency profile information, frequency domain sub-band coefficients under the j yardstick

Be expressed as successively:

$H_k^j=A_k^{j-1}*h_1^{j-1}*h_0^{j-1},$ $V_k^j=A_k^{j-1}*h_0^{j-1}*h_1^{j-1}$

$D_k^j=A_k^{j-1}*h_1^{j-1}*h_1^{j-1},$ $A_k^j=A_k^{j-1}*h_0^{j-1}*h_0^{j-1},$ j＝1，2，3

In the formula,

Be the low pass resolution filter coefficient of " 97 " small echo,

Be the high pass resolution filter coefficient of " 97 " small echo,

With

Be respectively right

With Carry out the filter coefficient that j interpolation obtains,

Represent subgraph to be transformed;

(1e) to high-frequency sub-band coefficient sequence W ^j=(H ^j, V ^j, D ^j) and low frequency sub-band coefficient sequence A ³Carrying out sampling matrix is M _ΛThe matrix interpolation computing, obtain the high-frequency sub-band coefficient sequence after the interpolation: W ^{J '}=(H ^{J '}, V ^{J '}, D ^{J '}) and low frequency sub-band coefficient sequence A ^{3 '}, be expressed as:

j＝1，2，3

In the formula, M _ΛBe sampling matrix, n=(n ₁n ₂) position coordinates of remarked pixel point in image,

Expression is carried out the new position coordinates that obtains after the matrix interpolation computing to the position coordinates of pixel, and Λ represents frequency domain sub-band coefficients W _kThe subscript collection;

(1f) respectively to the high-frequency sub-band coefficient sequence W behind the interpolation arithmetic ^{J '}=(H ^{J '}, V ^{J '}, D ^{J '}) and low frequency sub-band coefficient sequence A ^{3 '}Superpose, obtain input picture I _InCarry out the frequency domain representation coefficient D=(DH after the non-lower sampling Directionlet conversion ^j, DV ^j, DD ^j, DA ³),

Wherein: DH ^jBe meant the change direction high fdrequency component of detailed information of presentation video, DV ^jBe meant the high fdrequency component of presentation video formation direction detailed information, DD ^jBe meant that presentation video is changed direction and the high fdrequency component of formation direction detailed information, DA ³Be meant the low frequency component of presentation video low frequency profile information, they are expressed as respectively:

${\mathrm{DH}}^j\left(n\right)=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\&Lambda;}}\right)\right|-1}{\mathrm{\&Sigma;}}}{H}_k^{j^{\&prime;}}(n+S_k),$

${\mathrm{DV}}^j\left(n\right)=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\&Lambda;}}\right)\right|-1}{\mathrm{\&Sigma;}}}{V}_k^{j^{\&prime;}}(n+S_k),$

${\mathrm{DD}}^j\left(n\right)=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\&Lambda;}}\right)\right|-1}{\mathrm{\&Sigma;}}}{D}_k^{j^{\&prime;}}(n+S_k),$

${\mathrm{DA}}^3\left(n\right)=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\&Lambda;}}\right)\right|-1}{\mathrm{\&Sigma;}}}{A}_k^{3^{\&prime;}}(n+S_k),$ j＝1，2，3

In the formula, n=(n ₁n ₂) position coordinates of remarked pixel point in image, S _kRepresent k displacement vector, (n+S _k) expression new position coordinates that position coordinates is shifted and obtains, DH ^j, DV ^j, DD ^jForm the high fdrequency component among the frequency domain representation coefficient D

DA ³Form the low frequency component among the frequency domain representation coefficient D

3. according to claims 1 described mammary X-ray figure Enhancement Method, wherein the white Gaussian noise observing matrix Φ of the described use generation of step (3) is to the high fdrequency component among the frequency domain representation coefficient D

Observe, be meant observing matrix Φ that use generates at random and the high fdrequency component among the frequency domain representation coefficient D

Multiply each other, obtain observed reading

Wherein:

Be to high fdrequency component

In the change direction coefficient DH of detailed information of presentation video ^jThe result who observes is expressed as:

Be to high fdrequency component

The coefficient DV of middle presentation video formation direction detailed information ^jThe result who observes is expressed as:

Be to high fdrequency component

In presentation video change direction and the coefficient DD of formation direction detailed information ^jThe result who observes is expressed as:

4. according to claims 1 described mammary X-ray figure Enhancement Method, the described employing of step (4) OMP algorithm wherein, X recovers to observed reading, carries out as follows:

(4a) iteration count t is initialized as 1, i.e. t=1;

(4b) calculating observation matrix

In the column vector of each row

With observed reading

In i part coefficient

P row column vector v _pBetween inner product, obtain inner product sequence I, each inner product is I among the I _jBe expressed as:

p＝1，2，...，n，j＝1，2，...，N

In the formula, n is a coefficient Columns, N is the columns of observing matrix Φ;

(4c) inner product of numerical value maximum among the calculating inner product sequence I:

(4d) calculate maximum inner product I _MaxPosition number in inner product sequence I:

(4e) with maximum inner product I _MaxPosition number λ _tConstitute sequence number collection: Λ _t=(Λ _T-1, λ _t), in the formula, sequence number collection Λ ₀Be empty set:

The expression empty set;

(4f) with λ among the observing matrix Φ _tThe column vector of row

Constitute augmented matrix:

In the formula, augmented matrix Φ ₀Be empty set:

(4g) with column vector From observing matrix Φ, remove, be about to λ among the observing matrix Φ _tThe column vector of row

Be changed to empty set:

(4h) utilize augmented matrix Φ _tTo column vector v _pEstimate, obtain new estimated value x _t:

$x_t=\frac{<{\mathrm{\&Phi;}}_t,v_p>}{<{\mathrm{\&Phi;}}_t,{\mathrm{\&Phi;}}_t>}$

In the formula,＜Φ _t, v _pBe meant augmented matrix Φ _tWith column vector v _pInner product,＜Φ _t, Φ _tBe meant augmented matrix Φ _tWith the inner product of itself;

(4i) according to new estimated value x _tWith position number λ _tCalculate the column vector after reconstruct recovers:

(4j) utilize new estimated value x _t, augmented matrix Φ _tWith column vector v _pCalculate residual values r _t: r _t=v _p-Φ _tx _t

(4k) increase iteration count t:t=t+1, use residual values r _tReplacement step (4b) and (4h) in column vector v _p, repeating step (4b)～(4k) equals 1/4th of observation frequency M up to iteration count t, promptly

Till, the column vector that obtain this moment

P=1,2 ..., n has constituted through the high frequency coefficient after the reconstruct recovery Be expressed as

N is a coefficient

Columns, the high frequency coefficient after reconstruct recovers

Constituted the frequency domain high fdrequency component: This frequency domain high fdrequency component is the main high fdrequency component feature of representative image more concentratedly, and then improves the enhancing effect of subsequent treatment to image.

5. according to claims 1 described mammary X-ray figure Enhancement Method, wherein step (5) is described to the high fdrequency component after recovering

Carry out linear enhancement process, be meant the frequency domain high frequency coefficient after the reconstruct recovery

Multiply by amplification factor μ, the frequency domain high frequency coefficient after being enhanced is:

Be expressed as

The present invention is on the basis of carrying out repeatedly experiment test and checking, and the span that draws reinforcing coefficient μ is between 3～6.

6. according to claims 1 described mammary X-ray image enchancing method, wherein step (6) is described to the low frequency component among the frequency domain representation coefficient D

And the high fdrequency component after strengthening

Carry out non-lower sampling Directionlet inverse transformation, carry out as follows:

(6a) use sampling matrix M _ΛRespectively to the high fdrequency component after strengthening

With the low frequency component among the frequency domain representation coefficient D

Carry out the matrix sampling, obtain the frequency coefficient sequence

Wherein:

Be to use sampling matrix M _ΛTo high fdrequency component

Carry out the result of matrix sampling,

Be to use sampling matrix M _ΛTo low frequency component DA ³Carry out the result of matrix sampling, be expressed as

${\hat{Z}}_{\mathrm{ik}}^j\left(n\right)={\hat{Y}}_i^j\left(M_{\mathrm{\&Lambda;}}(n-S_k)\right),$ i＝1，2，3

$Z_k^3\left(n\right)={\mathrm{DA}}^3\left(M_{\mathrm{\&Lambda;}}(n-S_k)\right),$ k＝0，...，|det(M _Λ)|-1

In the formula, n=(n ₁n ₂) be the position coordinates of pixel in image, S _kRepresent k displacement vector, M _Λ(n-S _k) represent by displacement vector S _kPosition coordinates is shifted, and uses sampling matrix the position coordinates after being shifted to be carried out the new position coordinates that obtains after the matrix sampling;

(6b) with the frequency coefficient sequence

Carry out three layers of redundant wavelet reconstruction respectively, obtain the coefficient sequence Z behind the redundant wavelet reconstruction of j layer ^j, coefficient sequence Z ^jIn each subsystem number

Be expressed as:

$Z_k^j=Z_k^{j+1}*g_0^j*g_0^j+{\hat{Z}}_{3k}^{j+1}*g_1^j*g_1^j$

$+{\hat{Z}}_{2k}^{j+1}*g_1^j*g_0^j+{\hat{Z}}_{1k}^{j+1}*g_0^j*g_0^j,$ j＝2，1，0

In the formula,

With

Reconstruct low-pass filter coefficients and the reconstruct Hi-pass filter coefficient of representing " 97 " small echo respectively,

With

Be respectively right

With

Carry out the filter coefficient that j interpolation obtains; Obtain the image sequence after the conversion: Z this moment ⁰=Z ^j(j=0).

(6c) to the image sequence Z after the conversion ⁰In each subgraph

Carrying out sampling matrix is M _ΛThe matrix interpolation computing, obtain the image sequence Z after the interpolation, each subgraph Z among the image sequence Z _kBe expressed as:

In the formula, M _ΛBe sampling matrix, n=(n ₁n ₂) position coordinates of remarked pixel point in image,

The new position coordinates that obtains after the matrix interpolation is carried out in expression to position coordinates;

(6d) the image sequence Z after the matrix interpolation computing is superposeed, obtain through the enhancing image I after non-lower sampling Directionlet conversion and the compressed sensing enhancing _Out, be expressed as:

$I_{\mathrm{out}}=\underset{k=0}{\overset{\left|\det \left(M_{\mathrm{\&Lambda;}}\right)\right|-1}{\mathrm{\&Sigma;}}}{Z}_k(n+S_k)$

In the formula, n=(n ₁n ₂) position coordinates of remarked pixel point in image, S _kBe meant k displacement vector, Zk (n+S _k) represent image sequence Z _kIn k subgraph carry out the new subgraph that obtains after the coordinate displacement.

Images