CN102142133B - 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

Breast X-ray image enhancement method based on non-subsampled Directionlet transformation and compressed sensing

The technical field is as follows:

the invention belongs to the field of image processing, and relates to a non-subsampled Directionlet transform and compressed sensing image enhancement method.

Background art:

the breast cancer is one of the common malignant tumor diseases of women and seriously threatens the life health of the women. With the development of modern medical imaging technology, various new medical imaging technologies have been widely applied to various links such as medical diagnosis, preoperative planning, treatment, postoperative monitoring and the like, and these imaging technologies can comprehensively and accurately obtain various data of patients and provide more accurate information for diagnosis, treatment, operation and postoperative evaluation, wherein molybdenum target soft X-ray becomes the most reliable and commonly used means for diagnosing early breast cancer at present due to the advantages of higher spatial resolution, sensitivity to lumps and calcification, simple required equipment, low price and the like. However, the misdiagnosis and missed diagnosis rates when using breast molybdenum target soft X-ray films for diagnosis are still high, mainly due to poor image quality, benign manifestations of malignant lesions and visual fatigue or inattention of the observer. The inferior image quality is embodied in the following aspects: the contrast of the region of interest is poor, the intensity difference between the suspicious lesion region and the surrounding tissues is very weak, the shape of the lesion region is variable and has different sizes, the boundary is fuzzy, and the like. With the development of computers and related technologies, the problem can be effectively solved by enhancing the mammary gland X-ray image by using the computer technology, and doctors can be helped to understand and judge the image better.

In order to highlight the characteristics of the lesion area and improve the visual effect of the mammographic image, the most common method is to perform enhancement processing on the image. The traditional image enhancement processing technology achieves a good enhancement effect to a certain extent, but the enhancement effect of the traditional image enhancement processing technology cannot be satisfactory for mammary gland X-ray images with low contrast. To date, various enhancement methods have been proposed for mammographic images.

1. Unsharp masking method

The unsharp masking method is one of the commonly used enhancement methods in image processing, and is to add a certain proportion of high-frequency components of an image on the basis of an original image so as to achieve the effect of enhancing edge and detail information. The unsharp masking method has the advantages that the edge and detail information of the image can be well highlighted, and therefore the image enhancement effect is achieved. However, this unsharp masking method does not sufficiently consider the contrast information of the image, and therefore, when the contrast is low, the enhancement effect is not satisfactory.

2. Adaptive histogram equalization method

The adaptive histogram equalization method is a classical effective image enhancement method, and adopts a sliding window technology to perform histogram equalization on a window region containing a processed point, namely changing the histogram distribution of the region into uniform histogram distribution, and reassigns a pixel point to be processed in the window according to the mapping relation between the histogram and the gray level on the basis. The adaptive histogram equalization method has the advantages that the dynamic range of an image can be well adjusted, and meanwhile, the details of the image are enhanced, but the method also amplifies noise while improving the contrast, so that the enhancement effect of the method still needs to be improved.

3. Wavelet transformation adaptive gain processing method

The wavelet transform adaptive gain processing method is a relatively common image enhancement method, and includes the steps of firstly performing wavelet transform on an image, then performing adaptive enhancement processing according to the distribution condition of image transform coefficients, and finally performing inverse transform on the processed transform coefficients to obtain an enhanced image. The wavelet transform adaptive gain processing method has the advantages that the contrast of an image can be well improved, certain robustness is provided for noise, and when the gray level difference between a lesion region and a background region of the image is small, the enhancement effect of the image still needs to be improved.

Although the three methods have a good enhancement effect to a certain extent, the enhancement effect of the methods on the mammary X-ray image is not ideal because the mammary X-ray image has the characteristics of low contrast, rich noise, small gray scale difference between a lesion region and a background region and the like.

Disclosure of Invention

The invention aims to overcome the defects of the prior art and provides a mammary X-ray image enhancement method based on non-subsampled Directionlet transform and compressed sensing so as to effectively inhibit the background of an image, highlight the characteristics of a lesion area, improve the contrast of the image and make the enhancement effect of the mammary X-ray image more obvious.

The technical idea for realizing the aim of the invention is as follows: by removing the downsampling operation in the Directionlet transformation, the non-downsampling Directionlet transformation is realized, and the Directionlet transformation is combined with a compressive sensing method to enhance the clear effect on the mammary X-ray image. The specific scheme comprises the following steps:

(1) to input an image I_inConversion to the Directionlet transform domain, i.e. to the input image I using non-subsampled Directionlet transforms_inPerforming sub-band decomposition to obtain frequency domain representation coefficient D including high frequency components

And low frequency components

(2) And randomly generating a white gaussian noise observation matrix phi.

(3) Representing high frequency components in coefficients D to a frequency domain using a generated white Gaussian noise observation matrix phiAnd observing to obtain an observed value X.

(4) Recovering the observed value X by adopting an OMP algorithm to obtain a recovered high-frequency component

(5) For the recovered high frequency componentPerforming linear enhancement to obtain enhanced high-frequency component

(6) Representing the low frequency component in coefficient D for the frequency domain

And enhanced high frequency components

Performing non-subsampled Directionlet inverse transformation to obtain an enhanced image I_out。

The invention has the following advantages:

(1) the invention adopts the Directionlet conversion method, thereby effectively capturing the direction detail information of the image, obviously highlighting the detail characteristics of the lesion area in the mammary X-ray image and improving the definition of the image.

(2) The invention adopts the compressed sensing method to obtain the recovered high-frequency component, effectively concentrates the image energy, thereby intensively enhancing the high-frequency information and obviously improving the contrast of the image.

(3) The invention has simple structure and low calculation cost, and can simply and effectively enhance the image.

Description of the drawings:

FIG. 1 is a flow chart of the present invention.

FIG. 2 is a schematic diagram of the present invention for sampling an image using a sampling matrix.

FIG. 3 is a diagram of the results of the present invention after sampling an image using a sampling matrix.

FIG. 4 is a sub-flow diagram of the present invention for matrix interpolation and superposition of frequency domain coefficients using a sampling matrix.

FIG. 5 is a graph comparing the effect of the present invention on the treatment of a mammographic image with a conventional method.

The specific implementation scheme is as follows:

referring to fig. 1, the present invention includes non-downsampling Directionlet transform, compressed sensing, enhancement processing, and non-downsampling Directionlet inverse transform. The method comprises the following specific steps:

step 1: a sampling matrix is generated.

1.1) randomly generating two linearly independent integer vectors d starting from the origin of coordinates₁And d₂As the transformation direction and the queue direction of the image, respectively, wherein d₁＝[a₁，b₁]，d₂＝[a₂，b₂]，a₁Is an integer vector d₁Abscissa of the end point, b₁Is an integer vector d₁End point ordinate of a₂Is an integer vector d₂Abscissa of the end point, b₂Is an integer vector d₂The end point ordinate of (a);

1.2) integer vector d₁And d₂Form a sampling matrix M_Λ：

$M_{\mathrm{\&Lambda;}}=\left(\begin{array}{c}{d}_1\\ {\mathrm{<figure-callout\; id="2"\; label="d"\; filenames="HDA0000056203390000021.png"\; state="\{\{state\}\}">d</figure-callout>}}_2\end{array}\right)=\left(\begin{array}{cc}{a}_1& b_1\\ a_2& {\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="HDA0000056203390000021.png"\; state="\{\{state\}\}">b</figure-callout>}}_2\end{array}\right),a_1,a_2,b_1,{\mathrm{<figure-callout\; id="2"\; label="b"\; filenames="HDA0000056203390000021.png"\; state="\{\{state\}\}">b</figure-callout>}}_2\&Subset;Z$

Wherein Z is a whole set of integers.

Step 2: according to a sampling matrix M_ΛTo input an image I_inDivided into | det (M)_Λ) I each other irrelevant subgraph sequence F, each subgraph in F is F_k，k＝0，...，|det(M_Λ)|-1，

Wherein, | det (M)_Λ) Is the sampling matrix M_ΛDeterminant of (F)_kIs any subgraph, represented as: f_k(n)＝I_in(M_Λ(n-S_k))

Wherein n is (n)₁ n₂) Is the position coordinate of the pixel point in the image, F_k(n) represents a position coordinate of (n) in the kth sub-diagram₁ n₂) Pixel point of (2), S_kRepresenting the kth sub-graph F_kDisplacement vector of (1)_in(M_Λ(n-S_k) Is the kth sub-graph, which is based on the displacement vector S_kTo input an image I_inShifting the pixel points and then using the sampling matrix M_ΛThe sampling principle is shown in fig. 2, and integer vectors d are respectively selected from fig. 2₁＝[1，1]And d₁＝[-1，1]As the transform direction and the queue direction, the generated sampling matrix is:

using displacement vectors S_kFor input image I_inAfter the position coordinates are shifted, the position coordinates are shiftedUsing the generated sampling matrix M_ΛMultiplying the position coordinates of the shifted image to obtain a sampled sub-image sequence, as shown in FIG. 3, a sampling matrix M_ΛDividing the input image into two sub-images, wherein the black dots represent sub-image F₁White dot representation subgraph F₂。

And step 3: and (3) carrying out three-layer redundant wavelet decomposition on the sub-graph sequence F to obtain a high-frequency sub-band coefficient sequence of the sub-graph sequence F under each scale: w^j＝(H^j，V^j，D^j) J ═ 1,2,3 and low frequency subband coefficient sequence a³Sequence of high frequency subband coefficients W^jEach subband coefficient of:

k＝0，...，|det(M_Λ) I-1, low frequency subband coefficient sequence A³Each subband coefficient of:

k＝0，...，|det)M_Λ)|-1；

wherein,

means that a sub-graph F is represented after the j-th layer redundant wavelet decomposition is carried out_kThe frequency domain sub-band coefficients of the mid-horizontal direction detail information,

means that a sub-graph F is represented after the j-th layer redundant wavelet decomposition is carried out_kThe frequency domain sub-band coefficients of the vertical direction detail information,

means that a sub-graph F is represented after the j-th layer redundant wavelet decomposition is carried out_kThe frequency domain subband coefficients of the mid-diagonal detail information,

after the third layer of redundant wavelet decompositionRepresentation sub-diagram F_kFrequency domain sub-band coefficient of middle and low frequency general picture information, frequency domain sub-band coefficient under j scaleExpressed sequentially as:

$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,

are the low-pass decomposition filter coefficients of the "97" wavelet,

are the high pass decomposition filter coefficients of the "97" wavelet,

and

are respectively a pair

And

the filter coefficients obtained by the interpolation of j times are carried out,

representing the subgraph to be transformed.

And 4, step 4: according to a sampling matrix M_ΛFor high-frequency subband coefficient sequences W respectively^j＝(H^j，V^j，D^j) And a low frequency subband coefficient sequence A³Obtaining an input image I through matrix interpolation and superposition_inFrequency domain representation coefficient D after non-downsampling Directionlet transform (DH)^j，DV^j，DD^j，DA³)。

Referring to fig. 4, the specific implementation of this step is as follows:

4.1) matrix interpolation operation

For high-frequency subband coefficient sequences W respectively^j＝(H^j，V^j，D^j) And a low frequency subband coefficient sequence A³Making the sampling matrix M_ΛObtaining the interpolated high-frequency sub-band coefficient sequence by the matrix interpolation operation: w^j′＝(H^j′，V^j′，D^j′) And a low frequency subband coefficient sequence A^3′Expressed as:

j＝1，2，3

in the formula, M_ΛIs a sampling matrix, n ═ n (n)₁ n₂) Indicating the position coordinates of the pixel points in the image,

representing a new position coordinate obtained by performing matrix interpolation operation on the position coordinate of the pixel point, wherein Λ represents a subscript set of a frequency domain subband coefficient;

4.2) superposition of sub-band coefficient sequences

Respectively carrying out interpolation operation on the high-frequency sub-band coefficient sequences W^j′＝(H^j′，V^j′，D^j′) And a low frequency subband coefficient sequence A^3′Overlapping to obtain an input image I_inFrequency domain coefficient D ═ after non-downsampling Directionlet transform (DH)^j，DV^j，DD^j，DA³)，

Wherein: DH^jIs a high frequency component, DV, representing detail information of image transformation direction^jIs a high frequency component, DD, representing detail information in the direction of the image queue^jIs a high frequency component, DA, representing detailed information of image conversion direction and alignment direction³Refers to low frequency components representing low frequency overview information of an image, which are respectively represented as:

${\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;}}}{\mathrm{<figure-callout\; id="3"\; label="A\; k"\; filenames="HDA0000056203390000021.png"\; state="\{\{state\}\}">A</figure-callout>}}_k^{}(n+S_k),$ j＝1，2，3

wherein n is (n)₁ n₂) Representing the position coordinates of the pixel points in the image, S_kDenotes the kth displacement vector, (n + S)_k) Representing new position coordinates, DH, obtained by shifting the position coordinates^j，DV^j，DD^jHigh frequency components in the constituent frequency-domain representation coefficients D

DA³Composing the low-frequency component in the frequency-domain representation coefficient D

And 5: a Gaussian random matrix phi is randomly generated to serve as an observation matrix, the size of phi is M multiplied by N, M is the observation frequency, N is the number of rows of frequency domain representation coefficients needing to be observed, the value range of the observation frequency M is 250-500 on the basis of multiple experimental tests and verification, and the observation frequency of the example is 400.

Step 6: observing the high-frequency component in the frequency domain representation coefficient D by using the observation matrix phi, namely randomly generating the observation matrix phi and the high-frequency component in the frequency domain representation coefficient D

Multiplying to obtain observed value

Wherein:

is to the high frequency component

Coefficient DH representing detail information on image conversion direction^jThe results of the observations are expressed as:

is to the high frequency component

Coefficient DV representing detail information indicating image queue direction^jThe results of the observations are expressed as:

is to the high frequency componentCoefficient DD representing detail information of image transformation direction and queue direction^jThe results of the observations are expressed as:

and 7: restoring the observed value X by using the observation matrix phi generated in the step 2 and using an orthogonal matching pursuit OMP method, wherein the restored high-frequency coefficient is

7.1) initializing the iteration counter t to 1, i.e. t ═ 1;

7.2) calculating the Observation matrix

Column vector of each column in the column

And the observed value

Coefficient of the ith part

Of the p-th column vector v_pInner products between them to obtain inner product sequence I, each inner product in I is I_jExpressed as:

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

in which n is a coefficient

N is the number of columns of the observation matrix Φ;

7.3) calculating the largest numerical inner product in the inner product sequence I:

7.4) calculating the maximum inner product I_maxPosition number in inner product sequence I:

7.5) maximum inner product I_maxPosition number of (2)_tForm a sequence number set: lambda_t＝(Λ_t-1，λ_t) In the formula, the serial number set is Λ₀As an empty set:

representing an empty set;

7.6) using the lambda-th of the observation matrix phi_tColumn vector of a column

Constructing an augmentation matrix:in the formula, the augmentation matrix phi₀As an empty set:

7.7) vector the columns

Removing from the observation matrix phi, i.e. λ of the observation matrix phi_tColumn vector of a column

Setting as an empty set:

7.8) Using the augmentation matrix Φ_tVector v of alignment_pEstimating to obtain a new estimated value x_t：

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

In the formula,<Φ_t，v_p>is the finger augmentation matrix phi_tAnd column vector v_pThe inner product of (a) is,<Φ_t，Φ_t>is the finger augmentation matrix phi_tInner product with itself;

7.9) based on the new estimated value x_tAnd position number λ_tAnd (3) calculating a column vector after reconstruction recovery:

7.10) Using the new estimate x_tPhi, an amplification matrix_tAnd column vector v_pCalculating residual value r_t：r_t＝v_p-Φ_tx_t；

7.11) increment the iteration counter t: t is t +1, using residual r_tReplacement of the column vector v in steps 7.2) and 7.8)_pRepeating steps 7.2) to 7.11) until the iteration counter t is equal to one quarter of the number of observations M, i.e.

The column vector obtained at this time

N constitutes the high-frequency coefficient after reconstruction recovery

Is shown as

n is a coefficient

The number of columns of (a), reconstructing the recovered high-frequency coefficients

The frequency domain high frequency components are formed:

the frequency domain high-frequency components can represent the main high-frequency component characteristics of the image in a centralized manner, and the image enhancement effect of subsequent processing is further improved.

And 8: for the high frequency coefficient after reconstruction recovery

Multiplying by the amplification factor mu to obtain the enhanced high-frequency coefficient as:

is shown as

On the basis of carrying out multiple experimental tests and verifications, the value range of the obtained enhancement coefficient mu is 3-6, and the observation frequency of the example is 4.

And step 9: using a sampling matrix M_ΛRespectively for the enhanced high-frequency components

And low frequency components in the frequency domain representation coefficient D

Matrix sampling is carried out to obtain a frequency domain coefficient sequence

Wherein:

is to use a sampling matrix M_ΛFor high frequency component

As a result of the matrix sampling being performed,

is to use a sampling matrix M_ΛFor low frequency component DA³The result of the matrix sampling, 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

Wherein n is (n)₁ n₂) Is the position coordinate of the pixel point in the image, S_kDenotes the kth displacement vector, M_Λ(n-S_k) Representing the passing displacement vector S_kAnd shifting the position coordinates, and performing matrix sampling on the shifted position coordinates by using a sampling matrix to obtain new position coordinates.

Step 10: sequence of frequency domain coefficients

Respectively carrying out three-layer redundant wavelet reconstruction to obtain a coefficient sequence Z after the jth layer redundant wavelet reconstruction^jCoefficient sequence Z^jEach sub-coefficient of

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,

andrepresenting the reconstructed low-pass filter coefficients and the reconstructed high-pass filter coefficients of the "97" wavelet respectively,

and

are respectively a pair

And

carrying out j times of interpolation to obtain a filter coefficient; this results in a transformed image sequence: z⁰＝Z^j(j＝0)。

Step 11: according to a sampling matrix M_ΛFor transformed image sequences Z respectively⁰Each of which is a subgraph

Obtaining an enhanced image I after non-subsampled Directionlet transform and compressed sensing enhancement through matrix interpolation and superposition_out。

11.1) matrix interpolation

For image sequence Z⁰Each of which is a subgraphMaking the sampling matrix M_ΛObtaining an interpolated image sequence Z by matrix interpolation operation, wherein each subgraph Z in the image sequence Z_kExpressed as:

in the formula, M_ΛIs a sampling matrix, n ═ n (n)₁ n₂) Indicating the position coordinates of the pixel points in the image,representing a new position coordinate obtained by performing matrix interpolation on the position coordinate;

11.2) superposition of image sequences

Superposing the image sequence Z after the matrix interpolation operation to obtain an enhanced image I after the non-subsampled Directionlet transform and the compressed sensing enhancement_outExpressed as:

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

wherein n is (n)₁ n₂) Representing the position coordinates of the pixel points in the image, S_kIs the k-th displacement vector, Z_k(n+S_k) Representing a sequence of images Z_kThe k-th one ofAnd carrying out coordinate shift on the subgraph to obtain a new subgraph.

The advantages of the present invention can be further illustrated by the following simulation experiments:

1. simulation conditions

The test images used in the present invention are derived from mammograms in the MIAS database.

2. Emulated content

2.1) the present invention performs enhancement test on four sets of mammary X-ray images by using a sharpening mask method, an adaptive histogram equalization method, a wavelet transform adaptive gain processing method, and the present invention method to obtain four sets of enhanced images, as shown in fig. 5, wherein fig. 5(a) is a graph of four sets of original images from top to bottom, fig. 5(b) is a graph of enhancement results of fig. 5(a) using a conventional mask sharpening method from top to bottom, fig. 5(c) is a graph of enhancement results of fig. 5(a) using a conventional adaptive histogram equalization method from top to bottom, fig. 5(d) is a graph of enhancement results of fig. 5(a) using a conventional wavelet transform adaptive gain processing method from top to bottom, and fig. 5(e) is a graph of enhancement results of fig. 5(a) using the present invention from top to bottom. Comparing fig. 5(e) with fig. 5(b), fig. 5(c), and fig. 5(d), respectively, it can be seen that the present invention can effectively enhance the lesion region of the breast image, and better suppress the influence of the normal tissue belonging to the background in the image on the enhanced lesion region, so that the complex background region becomes smooth, thereby further increasing the information content and readability of the image.

2.2) the present invention detects the target to background contrast ratio TB based on variance_cTaking the value as a judgment basis, respectively testing four groups of mammary X-ray images enhanced by using a sharpening mask method, a self-adaptive histogram equalization method, a wavelet transformation self-adaptive gain processing method and the method of the invention to obtain four groups of contrast ratios TB_cThe values are shown in Table 1, where the first column is the four sets of original image sequences and the second column is TB after four sets of images were enhanced using a sharpening mask_cValue, third column isTB after four groups of images are enhanced by self-adaptive histogram equalization method_cValue, column four is for four groups of images after TB enhancement using wavelet transform adaptive gain processing_cThe fifth column is TB after four groups of images were enhanced using the method of the invention_cValue, wherein, the ratio TB_cThe values are expressed as:

TB_c＝δ_μ/σ

in the formula, delta_μThe difference between the average gray-scale ratio of the detected target and the background in the original image and the enhanced image, delta, is measured_μExpressed as:wherein

And

means for detecting the mean value of the target T and the background B,

and

refers to the mean of the enhanced image. Sigma measures the degree of reduction of the divergence of the gray levels of the detection target in the enhanced image relative to the detection target in the original image,in the formula

And

the variance of the detected object, and thus the detected object to background contrast ratio TB, in the original image and the enhanced image, respectively_cA larger value indicates a better enhancement of the image.

TABLE 1 enhancement results evaluation of TB_c

As can be seen from Table 1, TB of the invention after testing four sets of images_cThe value is obviously greater than that of the other three existing methods, so the enhancement effect of the method is superior to that of the three existing methods, namely a sharpening mask method, a self-adaptive histogram equalization method and a wavelet transformation self-adaptive gain processing method, in objective measurement.

In conclusion, the method and the device improve the contrast of the detected target and the background, and effectively enhance the detail information of the mammary gland image.

Claims (6)

1. A mammary X-ray image enhancement method based on non-subsampled Directionlet transform and compressed sensing comprises the following steps:

(1) to input an image I_inConversion to the Directionlet transform domain, i.e. to the input image I using non-subsampled Directionlet transforms_inPerforming sub-band decomposition to obtain frequency domain representation coefficient D including high frequency components

And low frequency components

(2) Randomly generating a Gaussian white noise observation matrix phi;

(3) representing high frequency components in coefficients D to a frequency domain using a generated white Gaussian noise observation matrix phi

Observing to obtain an observed value X;

(4) recovering the observed value X by adopting an OMP algorithm to obtain a recovered high-frequency component

(5) For the recovered high frequency component

Performing linear enhancement to obtain enhanced high-frequency component

(6) Representing the low frequency component in coefficient D for the frequency domain

And enhanced high frequency componentsPerforming non-subsampled Directionlet inverse transformation to obtain an enhanced image I_out。

2. The mammography X-ray image enhancement method of claim 1, wherein the non-downsampling Directionlet transform is used to input image I in step (1)_inAnd (3) carrying out sub-band decomposition according to the following steps:

(1a) randomly generating two linearly independent integer vectors d₁And d₂Is divided intoRespectively as the transformation direction and the queue direction of the image, wherein d₁＝[a₁,b₁]，d₂＝[a₂,b₂]，a₁Is an integer vector d₁Abscissa of the end point, b₁Is an integer vector d₁End point ordinate of a₂Is an integer vector d₂Abscissa of the end point, b₂Is an integer vector d₂The end point ordinate of (a);

(1b) the integer vector d₁And d₂Form a 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$

Wherein Z is a full set of integers;

(1c) according to a sampling matrix M_ΛTo input an image I_inDivided into | det (M)_Λ) I each other irrelevant subgraph sequence F, each subgraph in F is F_k,k＝0,..,|det(M_Λ)|-1，

Wherein, | det (M)_Λ) Is the sampling matrix M_ΛDeterminant of (F)_kIs any subgraph, represented as:

F_k(n)＝I_in(M_Λ(n-S_k))

wherein n is (n)₁ n₂) Is the position coordinate of the pixel point in the image, F_k(n) represents a position coordinate of (n) in the kth sub-diagram₁ n₂) Pixel point of (2), S_kRepresenting the kth sub-graph F_kDisplacement vector of (1)_in(M_Λ(n-S_k) Is the kth sub-graph, which is based on the displacement vector S_kTo input an image I_inShifting the pixel points and then using the sampling matrix M_ΛThe matrix is obtained after matrix sampling operation is carried out on the matrix;

(1d) and (3) carrying out three-layer redundant wavelet decomposition on the sub-graph sequence F to obtain a high-frequency sub-band coefficient sequence of the sub-graph sequence F under each scale: w^j＝(H^j,V^j,D^j) J ═ 1,2,3 and low frequency subband coefficient sequence a³Sequence of high frequency subband coefficients W^jEach subband coefficient of: $W_k^j=(H_k^j,V_k^j,D_k^j),k=0,...,\left|\det \left(M_{\mathrm{\&Lambda;}}\right)\right|-1,$ low frequency subband coefficient sequence A³Each subband coefficient of:

wherein,

means that a sub-graph F is represented after the j-th layer redundant wavelet decomposition is carried out_kThe frequency domain sub-band coefficients of the mid-horizontal direction detail information,

means that a sub-graph F is represented after the j-th layer redundant wavelet decomposition is carried out_kThe frequency domain sub-band coefficients of the vertical direction detail information,

means that a sub-graph F is represented after the j-th layer redundant wavelet decomposition is carried out_kThe frequency domain subband coefficients of the mid-diagonal detail information,means that sub-graph F is represented after the third layer redundant wavelet decomposition_kFrequency domain sub-band coefficient of middle and low frequency general picture information, frequency domain sub-band coefficient under j scale

Expressed sequentially as:

$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=\mathrm{1,2,3}$

in the formula,

are the low-pass decomposition filter coefficients of the "97" wavelet,

are the high pass decomposition filter coefficients of the "97" wavelet, $h_0^{j-1},j=\mathrm{2,3}$ and $h_1^{j-1},j=\mathrm{2,3}$ are respectively a pair $h_0^{j-1},j=1$ And $h_1^{j-1},j=1$ the filter coefficient obtained by j-1 times of interpolation,representing a subgraph to be transformed;

(1e) for high frequency subband coefficient sequence W^j＝(H^j,V^j,D^j) And a low frequency subband coefficient sequence A³Making the sampling matrix M_ΛObtaining the interpolated high-frequency sub-band coefficient sequence by the matrix interpolation operation: w^j′＝(H^j′,V^j′,D^j′) And a low frequency subband coefficient sequence A^3′Expressed as:

in the formula, M_ΛIs a sampling matrix, n ═ n (n)₁ n₂) Indicating the position coordinates of the pixel points in the image,

representing new position coordinates obtained by performing matrix interpolation operation on the position coordinates of the pixel points, and Λ representing frequency domain sub-band coefficient W_kA subscript set of (a);

(1f) respectively carrying out interpolation operation on the high-frequency sub-band coefficient sequences W^j′＝(H^j′,V^j′,D^j′) And a low frequency subband coefficient sequence A^3′Overlapping to obtain an input image I_inFrequency domain representation coefficient D after non-downsampling Directionlet transform (DH)^j,DV^j,DD^j,DA³)，

Wherein: DH^jFinger watchHigh-frequency components, DV, of detail information of image transformation direction^jIs a high frequency component, DD, representing detail information in the direction of the image queue^jIs a high frequency component, DA, representing detailed information of image conversion direction and alignment direction³Refers to low frequency components representing low frequency overview information of an image, which are respectively represented as:

${\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;}}}{D}_k^{3\&prime;}(n+S_k),j=\mathrm{1,2,3}$

wherein n is (n)₁ n₂) Representing the position coordinates of the pixel points in the image, S_kDenotes the kth displacement vector, (n + S)_k) Representing new position coordinates, DH, obtained by shifting the position coordinates^j,DV^j,DD^jHigh frequency components in the constituent frequency-domain representation coefficients D

DA³Composing the low-frequency component in the frequency-domain representation coefficient D $\stackrel{\&OverBar;}{D}={\mathrm{DA}}^3.$

3. The mammary X-ray image enhancement method according to claim 1, wherein the step (3) of representing the high frequency components in the coefficients D to the frequency domain using the generated white Gaussian noise observation matrix Φ

Observing means using a randomly generated observation matrix phi and a high-frequency component in the frequency domain representation coefficient DMultiplying to obtain observed value $X=(X_1^j,X_2^j,X_3^j),$ Wherein:

is to the high frequency component

Coefficient DH representing detail information on image conversion direction^jThe results of the observations are expressed as: $X_1^j=\mathrm{\&Phi;}\&times;{\mathrm{DH}}^j;$

is to the high frequency component

Coefficient DV representing detail information indicating image queue direction^jThe results of the observations are expressed as: $X_2^j=\mathrm{\&Phi;}\&times;{\mathrm{DV}}^j;$

is to the high frequency component

Coefficient DD representing detail information of image transformation direction and queue direction^jThe results of the observations are expressed as:

4. the mammary X-ray image enhancement method according to claim 1, wherein the step (4) of recovering the observed value X by using an OMP algorithm is performed by the following steps:

(4a) initializing an iteration counter t to 1, namely t is 1;

(4b) calculating an observation matrix

Column vector of each column in the column

And the observed value

Coefficient of the ith part

Of the p-th column vector v_pInner products between them to obtain inner product sequence I, each inner product in I is I_jExpressed as:

in which n is a coefficient

N is the number of columns of the observation matrix Φ;

(4c) computingInner product sequence I inner product with the largest value:

(4d) calculating the maximum inner product I_maxPosition number in inner product sequence I:

(4e) by maximum inner product I_maxPosition number of (2)_tForm a sequence number set: lambda_t＝(Λ_t-1,λ_t) In the formula, the serial number set is Λ₀As an empty set: representing an empty set;

(4f) by lambda of the observation matrix phi_tColumn vector of a column

Constructing an augmentation matrix:

in the formula, the augmentation matrix phi₀As an empty set:

(4g) vector the column

Removing from the observation matrix phi, i.e. λ of the observation matrix phi_tColumn vector of a column

Setting as an empty set:

(4h) using an amplification matrix phi_tVector v of alignment_pEstimating to obtain a new estimated value x_t：

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

In the formula,<Φ_t,v_p>is the finger augmentation matrix phi_tAnd column vector v_pThe inner product of (a) is,<Φ_t,Φ_t>is the finger augmentation matrix phi_tInner product with itself;

(4i) based on the new estimated value x_tAnd position number λ_tAnd (3) calculating a column vector after reconstruction recovery:

(4j) using the new estimate x_tPhi, an amplification matrix_tAnd column vector v_pCalculating residual value r_t：r_t＝v_p-Φ_tx_t；

(4k) Increment an iteration counter t: t is t +1, using residual r_tReplacing the column vector v in steps (4b) and (4h)_pRepeating steps (4b) - (4k) until the iteration counter t is equal to one fourth of the number of observations M, i.e.

The column vector obtained at this time

Form the high-frequency coefficient after reconstruction and recovery

Is shown as

n is a coefficient

The number of columns of (a), reconstructing the recovered high-frequency coefficients

The frequency domain high frequency components are formed:

the frequency domain high-frequency components can represent the main high-frequency component characteristics of the image in a centralized manner, and the image enhancement effect of subsequent processing is further improved.

5. The mammography X-ray image enhancement method of claim 1, wherein the pair of the restored high-frequency components of step (5)

Performing linear enhancement processing on the reconstructed frequency domain high-frequency coefficient

Multiplying by an amplification factor mu to obtain an enhanced frequency domain high-frequency coefficient as follows:

is shown asThe value range of the amplification factor mu is 3-6.

6. The mammary X-ray image enhancement method according to claim 1, wherein the step (6) of representing the low frequency component in the coefficient D to the frequency domain

And enhanced high frequency components

Performing non-downsampling Directionlet inverse transformation according to the following steps:

(6a) using a sampling matrix M_ΛRespectively for the enhanced high-frequency components

And low frequency components in the frequency domain representation coefficient DMatrix sampling is carried out to obtain a frequency domain coefficient sequence

Wherein:

is to use a sampling matrix M_ΛFor high frequency component

As a result of the matrix sampling being performed,

is to use a sampling matrix M_ΛFor low frequency component DA³The result of the matrix sampling, expressed as:

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

$Z_k^3\left(n\right)=D{A}^3\left(M_{\mathrm{\&Lambda;}}(n-S_k)\right),k=0,...,\left|\det \left(M_{\mathrm{\&Lambda;}}\right)\right|-1$

wherein n is (n)₁ n₂) Is the position coordinate of the pixel point in the image, S_kRepresents the kth bitMotion vector, M_Λ(n-S_k) Representing the passing displacement vector S_kShifting the position coordinates, and performing matrix sampling on the shifted position coordinates by using a sampling matrix to obtain new position coordinates;

(6b) sequence of frequency domain coefficients

Respectively carrying out three-layer redundant wavelet reconstruction to obtain a coefficient sequence Z after the jth layer redundant wavelet reconstruction^jCoefficient sequence Z^jEach sub-coefficient of

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=\mathrm{2,1,0}$

in the formula,

and

representing the reconstructed low-pass filter coefficients and the reconstructed high-pass filter coefficients of the "97" wavelet respectively, $g_0^j(j=\mathrm{1,2})$ and $g_1^j(j=\mathrm{1,2})$ are respectively a pair $g_0^j(j=0)$ And $g_1^j(j=0)$ carrying out j times of interpolation to obtain a filter coefficient; this results in a transformed image sequence: z⁰＝Z^j(j＝0)；

(6c) For the transformed image sequence Z⁰Each of which is a subgraph

Making the sampling matrix M_ΛObtaining an interpolated image sequence Z by matrix interpolation operation, wherein each subgraph Z in the image sequence Z_kExpressed as:

in the formula, M_ΛIs a sampling matrix, n ═ n (n)₁ n₂) Indicating the position coordinates of the pixel points in the image,representing a new position coordinate obtained by performing matrix interpolation on the position coordinate;

(6d) superposing the image sequence Z after the matrix interpolation operation to obtain an enhanced image I after the non-subsampled Directionlet transform and the compressed sensing enhancement_outExpressed as:

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

wherein n is (n)₁ n₂) Representing the position coordinates of the pixel points in the image, S_kIs the k-th displacement vector, Z_k(n+S_k) Representing a sequence of images Z_kAnd carrying out coordinate shift on the kth sub-graph to obtain a new sub-graph.

Images