CN103985109B - Feature-level medical image fusion method based on 3D (three dimension) shearlet transform - Google Patents

Abstract

The invention discloses a feature-level medical image fusion method based on 3D (three dimension) shearlet transform, belonging to the technical field of medical image processing and application. The feature-level medical image fusion method mainly comprises the following steps: 1, performing 3D-D-CSST (three dimensional-discrete-compact shearlet transform) or 3D-DT-CSST (three dimensional dual-tree compact shearlet transform) on two images to obtain transformation coefficient images Ca and Cb; 2, performing image fusion on transformation coefficients to obtain a fusion coefficient Cf; and 3, performing DWT or DTCWT inverse transformation, performing backward shear transformation on the transformed image to obtain a fusion image Vf. According to the feature-level medical image fusion method, the problems that the quality of a fused image is relatively low and information which is partially important and is not remarkable is easily ignored are solved.

Description

A kind of feature level Method of Medical Image Fusion that wave conversion is sheared based on 3D

Technical field

The invention belongs to Medical Image Processing and applied technical field, and in particular to a kind of spy that wave conversion is sheared based on 3D Levy a grade Method of Medical Image Fusion, solve that fused image quality is relatively low and local is important but inapparent information is easily neglected The problem omitted.

Background technology

Medical image fusion is one kind of image co-registration, and many methods have been widely used in clinical diagnosis.Fusion Referring to will become a width as the important information in source images that the distinct devices such as CT, MRI are gathered with regard to target is extracted and merged The process of image.The information included in the generated image of distinct device or the different configurations of same equipment is different, some letters Breath has similitude, but most information is complementary.For example, the mainly human body that CT images are provided is dense, hard tissues Information, and MRI image then mainly provides the information of soft tissue.The post processing figure of the same information for once gathering of same MRI machine Picture, such as T2* provide the comparative information of tissue relaxation time, quantity of magnetism figure (QSM：Quantitative Susceptibility Mapping the magnetic susceptibility comparative information caused by various magnetic bio labels (such as iron, calcium, gadolinium contrast medium)) is provided. Typically, image co-registration is needed first by source image registration, and T2* and QSM images are to be based on to be carried out with the data of single pass Post processing is generated, so the two is completely registering.

The research of current medical image co-registration primary concern is that the situation of two dimensional image, but multiclass Medical Devices now All it is to generate 3-D view.In 3-D view each point gray value not only with same layer neighbor point cross-correlation, also with adjacent layer in Neighbor point cross-correlation.Traditional two-dimentional fusion method can cause the loss of third dimension information, it is therefore necessary to which research can be directly Process the fusion method of 3-D view.

Blending algorithm can be processed in spatial domain or transform domain.In spatial domain, fused images are typically the weighting of source data Averagely, such method is simply easily achieved, but fused image quality is not high.Transform domain method follows following steps：1) source figure As transforming to transform domain, 2) image coefficient is processed by fusion criterion, the coefficient after being merged, 3) finally become coefficient Spatial domain is returned to, output is fused images.Research emphasis are concentrated mainly at 2 points in this kind of algorithm：The selection and fusion of conversion is accurate Design then.Many multiple dimensioned (Multi-Scale) conversion can be applied in blending algorithm, such as DWT, DTCWT, Curvelet, Shearlet etc..

Shearing wave conversion is to propose in recent years and the efficient conversion for representing of progressively ripe multidimensional data.In fact, being directed to Wavelet transformation lacks shortcoming to edge isotropy feature rarefaction representation, and scholars are also suggested many other multiple dimensioned Conversion.But it is unique while possessing the conversion of advantages below in all methods to shear wave conversion：Only one of which or limited generation Function set, expression high dimensional data that can be almost optimum is uniformly processed to continuous data and discrete data, possesses tight Zhi Shixian etc. Deng.Shearing wave conversion is widely used in image procossing, such as denoising, edge detection, strengthen etc..

Shearing wave is applied equally to image co-registration, and existing image fusion technology has the disadvantage that：1st, tradition is based on The fusion method of wavelet transformation and pyramid transform, because multi-scale transform lacks the rarefaction representation energy to picture structure directionality Power, causes the quality of fused images relatively low；2nd, the image co-registration based on Pixel-level, does not account for the structural information of image, When can cause image co-registration, the important really inapparent information in local and be ignored.These defects can be to final medical diagnosis Have a negative impact.

The content of the invention

For above-mentioned prior art, present invention aim at providing a kind of feature level medical science figure that wave conversion is sheared based on 3D As fusion method, the fusion method based on wavelet transformation and pyramid transform is solved, because multi-scale transform lacks to picture structure The rarefaction representation ability of directionality and cause the quality of fused images relatively low；And during image co-registration, local is important but not The defect that significant information is easily ignored, these defects eventually have a negative impact to medical diagnosis.

In order to solve above-mentioned technical problem, the present invention is adopted the following technical scheme that：

Herein 3D shearing waves specifically refer to 3D and tightly prop up tight shearing wave (3D- of shearing wave (3D-D- shearing waves) or the double trees of 3D DT- shearing waves), D- shearing wave conversions include two steps：Forward direction shear is converted and DWT conversion；DT- shearing wave conversions are included Two steps：Forward direction shear is converted and DTCWT conversion.

A kind of feature level Method of Medical Image Fusion that wave conversion is sheared based on 3D, it is characterised in that comprise the steps：

First, two width 3D medical image V to be fused are prepared_a、V_b, three directions of two width images are carried out respectively before to Shear is converted, and wavelet transform DWT or bi-input bi-output system conversion DTCWT is carried out to the image after conversion, obtains corresponding Multigroup changing image coefficient C_a、C_b；

2nd, image co-registration is carried out to the coefficient that 3D shearing wave conversions are obtained, obtains fused images coefficient C_f；

3rd, to image coefficient C after step 2 fusion_fDWT or DTCWT inverse transformations are carried out, the image after conversion is entered The backward shear conversion of row obtains multigroup fused images, and to these images final fused images V are averagely obtained_f。

In the present invention, the detailed step of the step 2 includes following two step：

2.1st, the image C that wave conversion is obtained is sheared to 3D_a、C_bLow frequency part C_aL、C_bLMerged using average criterion Low frequency part C of image_fL；

2.2nd, to HFS C_aH、C_bHUsing the fusion of feature level, the feature class of same position image to be fused is judged Type, is merged by the maximum information criterion that retains, and obtains C_fH；

In the present invention, in the step one, to shear conversion before first carrying out to image, then to the image after conversion Carry out wavelet transform DWT or bi-input bi-output system conversion DTCWT；Forward direction shear conversion is specific as follows：For one group three Dimension data l × m × n sets up coordinate system, and origin is (0,0,0), and its angle steel joint is (l-1, m-1, n-1), and three sides are carried out to it To shear conversion it is as follows：Shear conversion wherein for z directions is referred to and carries out following coordinate change to the point in data Change：

It is for the shear transformation for mula in x directions：

Doing shear transformation for mula for y directions is：

Wherein, (x, y, z) is the coordinate before conversion, and (x ', y ', z ') is the coordinate after conversion.k_tr, tr=a1, b1, a2, B2, a3, b3 } it is mobile ultimate range.k_trDifferent values are taken, the information for retaining different directions, thus shearing wave will be obtained Conversion can be producedIndividual 3D rendering, whereinWithFor k_aiAnd k_biDirection number.

Feature-based fusion is adopted to the low frequency part of changing image in the step 2.1, its fusion rule is：

Cf_L=(C_aL+C_aL)/2(2)

Feature-based fusion is adopted to the HFS of changing image in the step 2.2, concrete operation step is as follows：

2.2.1, changing image coefficient C is first calculated_a、C_bHFS C_aH、C_bHStructure tensor, then structure tensor is entered Row rank is analyzed：

For changing image coefficient C_a、C_bHFS C_aH、C_bHEach point, structure tensor is 3 × 3 matrixes, Rank of matrix desirable 0,1,2,3, flat, planar respectively in correspondence image, wire, dotted region feature；Ω is regional area l₁×m₁×n₁, the structure tensor of point p is expressed as

W (r) is a l₁×m₁×n₁The Gaussian template of size；V_x(p)、V_y(p)、V_zP () is respectively image to x, y, z axle Partial derivative on three directions；

Calculate characteristic value E of this 3 × 3 tensor matrix_x、E_y、E_z, given threshold K is control parameter, is set to 0.01, the nonzero eigenvalue number of point pFor the same position of two width figures, C is remembered_aNonzero eigenvalue number be M_a, Note C_bNonzero eigenvalue number be M_b, M_a、M_bAs the approximate of tensor rank of matrix；

If 2.2.2, M_a=M_b, then two width figures have same type feature in this position, calculate the phase of this position Like degree

Calculate threshold valueFusion rule is：

γ_abDuring≤α, this position is redundancy, selects weighted criterion：

C_fH=ω_aC_aH+ω_bC_bH (5)

γ_ab>During α, this position is complementary information, using MRE criterions:

If 2.2.3, M_a≠M_b, fusion criterion：

The step 3 does backward shear conversion to the image after DWT or DTCWT inverse transformations, specific as follows：

To the inverse operation of shear conversion before referring to for backward shear conversion, wherein the shear conversion for z directions is Finger carries out following coordinate transform to the point in data：

It is for the shear transformation for mula in x directions：

Doing shear transformation for mula for y directions is：

Wherein, (x, y, z) be conversion before coordinate, (x ', y ', z ') be conversion after coordinate, k_tr, tr=a1, b1, a2, B2, a3, b3 } value corresponding to front to shear conversion institute values.

Compared with prior art, the invention has the advantages that：

First, relative to the rarefaction representation energy lacked based on tradition Wavelet and pyramid transform etc. to direction architectural characteristic For the multi-scale transform of power, compact schemes shearing wave conversion has the almost optimum each energy to different feature represented in high dimensional signal Power, fused images retain more accurately directional information, cause fusion mass higher；

2nd, DT- shearings wave conversion introduces double tree constructions, reduces and moves the fused images distortion that denaturation is caused；

3rd, the spatial domain compact sup-port of DT- shearing waves and D- shearing waves, relative to frequency domain shearing wave, fusion mass is higher；

4th, the present invention uses feature level fusing method, it is contemplated that scanning organ internal structural characteristic (including flat, planar, Wire, dotted region), to greatest extent object of reservation structural information and physical features, special relative to only consideration high frequency coefficient statistics The Pixel-level fusion criterion levied, fused image quality is higher；

5th, the quality of fused images of the present invention is weighed by objective indicator (MI, QAB | F), and its quality is higher.

Description of the drawings

Fig. 1 is image interfusion method schematic diagram of the present invention；

Fig. 2 is that two dimension shear converts schematic diagram；

Fig. 3 is that three-dimensional shear converts (z-axis direction shear conversion) schematic diagram.

Specific embodiment

Below in conjunction with the drawings and the specific embodiments, the invention will be further described.

By taking T2* magnitude images and QSM images as an example, this experiment is schemed with the inventive method to three-dimensional T2* magnitude images and QSM As being processed, fused images are finally given, in example, image size is 128 × 128 × 128.

The image interfusion method of the present embodiment considers how first the anisotropy for showing 3-D view, and shear conversion can With each to different feature of good behaviour image；Secondly consider the impact of the translation qualitative change that DWT conversion brings, thus carry out double Tree Complex Wavelet Transform DTCWT；Finally consider to be fused to the information of low frequency coefficient and high frequency coefficient institute band is as much as possible In closing coefficient image：For low frequency adopts average criterion；For high frequency coefficient is merged using the fusion criterion of feature level.

Flow process is as shown in figure 1, comprise the following steps：

Step one：To two width image V_a、V_bConversion coefficient C is obtained before carrying out to 3D-DT- shearing wave conversions_a、C_b.Performed Journey is included before three-dimensional to shear conversion and three-dimensional DTCWT conversion.

Forward direction shear is converted：For one group of three-dimensional data l × m × n sets up coordinate system, origin is (0,0,0), and its is diagonal Point is (l-1, m-1, n-1), shear conversion in z-axis direction is carried out to it and refers to that x to the point in data, y-coordinate enter line translation

It is to the shear transformation for mula of other both directions：

Fig. 2 is the schematic diagram of shear conversion, for this example, l=m=n=128, selects k_tr, tr=a1, b1, a2, B2, a3, b3 }=- 64,0,64,27 groups of changing image data will be produced；

Step 2：27 system number view data C are obtained to 3D-DT- shearing wave conversions_a、C_bCarry out fusion and obtain C_f.Perform Process includes merging low frequency coefficient and merging high frequency coefficient.

The coefficient C that wave conversion is obtained is sheared to 3D-DT-_a、C_bLow frequency part C_aL、C_bLMerged using average criterion：

C_fL=(C_aL+C_bL)/2

2) for HFS C_aH、C_bH, the regional area Ω sizes of selected element p are l₁×m₁×n₁, elect 3 × 3 as here × 3, calculate the structure tensor of p points：

W (r) is a l₁×m₁×n₁The Gaussian template of size；V_x(p)、V_y(p)、V_zP () is respectively image to x, y, z axle Partial derivative on three directions.

There is correlation between voxel, it is difficult to have proper flat, planar, wire and dotted region, so When extracting image spatial feature, to nonzero eigenvalue, this condition has made appropriate relaxing.For a certain characteristic value of structure tensor During less than respective threshold, it is believed that this characteristic value is zero, and then think that more than the characteristic value number of threshold value be rank of matrix.Calculate this Characteristic value E of 3 × 3 matrixes_x、E_y、E_z, given thresholdk For control parameter, 0.01 can be set to, the nonzero eigenvalue number of point pM is approximately and opens The order of moment matrix, for the same position of two width figures, M_aRecord C_aHMore than the characteristic value number of threshold values, M_bRecord C_bHMore than threshold values Characteristic value number；

If M_a=M_b, then two width figures have same type feature in this position, then calculate the similarity of this position

Calculate threshold valueFusion rule is：

γ_abDuring≤α, this position is redundancy, selects weighted criterion

C_fH=ω_aC_aH+ω_bC_bH

γ_ab>During α, this position is complementary information, using MRE criterions:

If M_a≠M_b, fusion criterion：

Step 3：To fusion coefficients image C_fCarry out backward 3D-DT- shearings wave conversion and obtain final fused images.Perform Process is included before three-dimensional to shear conversion and three-dimensional DTCWT conversion

To the inverse operation of shear conversion before referring to for backward shear conversion, wherein the shear conversion for z directions is Finger carries out following coordinate transform to the point in data：

It is for the shear transformation for mula in x directions：

Doing shear transformation for mula for y directions is：

L=m=n=128 is selected, k is selected_tr, { tr=a1, b1, a2, b2, a3, b3 }=- 64,0,64；Finally to 27 3D rendering after reciprocal transformation does and averagely obtain fused images V_f。

The above, the only preferred embodiment of the invention, but protection scope of the present invention is not limited thereto, any ripe Those skilled in the art are known in scope disclosed in this invention, technology according to the present invention scheme and its inventive concept Equivalent or change in addition, belongs to protection scope of the present invention.

Claims (5)

1. it is a kind of based on 3D shear wave conversion feature level Method of Medical Image Fusion, it is characterised in that comprise the steps：

First, two width 3D medical image V to be fused are prepared_a、V_b, become to shear before carrying out to three directions of two width images respectively Change, wavelet transform DWT or bi-input bi-output system conversion DTCWT is carried out to the image after conversion, obtain corresponding multigroup Changing image coefficient C_a、C_b；

2nd, image co-registration is carried out to the coefficient that 3D shearing wave conversions are obtained, obtains fused images coefficient C_f；Concretely comprise the following steps：

2.1st, image coefficient C that wave conversion is obtained is sheared to 3D_a、C_bLow frequency part C_aL、C_bLMerged using average criterion Low frequency part C of image_fL；

2.2nd, image coefficient C that wave conversion is obtained is sheared to 3D_a、C_bHFS C_aH、C_bHMerged, obtained fused images HFS C_fH；

2.3rd, according to low frequency part C_fLWith HFS C_fHObtain fused images coefficient C_f；

3rd, to image coefficient C after step 2 fusion_fDWT or DTCWT inverse transformations are carried out, the image after conversion is carried out backward Shear conversion obtains fused images, and to these images final fused images V are averagely obtained_f。

2. it is according to claim 1 based on 3D shear wave conversion feature level Method of Medical Image Fusion, it is characterised in that In the step one, to shear conversion before first carrying out to image, then wavelet transform DWT is carried out to the image after conversion Or bi-input bi-output system conversion DTCWT；Forward direction shear conversion is specific as follows：Sit for one group of three-dimensional data l × m × n sets up Mark system, origin is (0,0,0), and its angle steel joint is (l-1, m-1, n-1), and the shear that three directions are carried out to it converts following institute Show：Shear conversion wherein for z directions is referred to and carries out following coordinate transform to the point in data：

It is for the shear transformation for mula in x directions：

Doing shear transformation for mula for y directions is：

Wherein, (x, y, z) is the coordinate before conversion, and (x ', y ', z ') is the coordinate after conversion；k_tr, tr=a1, b1, a2, b2, A3, b3 } it is mobile ultimate range；k_trDifferent values are taken, the information for retaining different directions will be obtained, thus shear wave conversion Can produceIndividual 3D rendering, whereinWithFor k_aiAnd k_biDirection number.

3. it is according to claim 1 based on 3D shear wave conversion feature level Method of Medical Image Fusion, it is characterised in that Two width 3D rendering V in the step 2.2_a、V_bCorresponding changing image coefficient C_a、C_bHFS C_aH、C_bHFusion method For：Using the fusion of feature level, the characteristic type of same position image to be fused is judged, carried out by the maximum information criterion that retains Fusion, obtains C_fH。

4. it is according to claim 3 based on 3D shear wave conversion feature level Method of Medical Image Fusion, it is characterised in that The HFS to changing image adopts feature-based fusion, and concrete operation step is as follows：

2.2.1, changing image coefficient C is first calculated_a、C_bHFS C_aH、C_bHStructure tensor, then row rank is entered to structure tensor Analysis：

For changing image coefficient C_a、C_bHFS C_aH、C_bHEach point, structure tensor is 3 × 3 matrixes, matrix Order desirable 0,1,2,3, flat, planar respectively in correspondence image, wire, dotted region feature；Ω is regional area l₁×m₁ ×n₁, the structure tensor of point p is expressed as

$\begin{array}{c}s(x,y,z)=\\ \left[\begin{array}{ccc}\underset{r\∈\mathrm{\Ω}}{\Σ}w\left(r\right)V_x^2(p-r)& \underset{r\∈\mathrm{\Ω}}{\Σ}w\left(r\right)V_x(p-r)V_y(p-r)& \underset{r\∈\mathrm{\Ω}}{\Σ}w\left(r\right)V_x(p-r)V_z(p-r)\\ \underset{r\∈\mathrm{\Ω}}{\Σ}w\left(r\right)V_x(p-r)V_y(p-r)& \underset{r\∈\mathrm{\Ω}}{\Σ}w\left(r\right)V_y^2(p-r)& \underset{r\∈\mathrm{\Ω}}{\Σ}w\left(r\right)V_y(p-r)V_z(p-r)\\ \underset{r\∈\mathrm{\Ω}}{\Σ}w\left(r\right)V_x(p-r)V_z(p-r)& \underset{r\∈\mathrm{\Ω}}{\Σ}w\left(r\right)V_y(p-r)V_z(p-r)& \underset{r\∈\mathrm{\Ω}}{\Σ}w\left(r\right)V_z^2(p-r)\end{array}\right]\end{array}---\left(2\right)$

W (r) is a l₁×m₁×n₁The Gaussian template of size；V_x(p)、V_y(p)、V_zP () is respectively image to three, x, y, z axle Partial derivative on direction；

Calculate characteristic value E of this 3 × 3 tensor matrix_x、E_y、E_z, given threshold K is control parameter, is set to 0.01, nonzero eigenvalue number M=su of point p_tm(E_t＞ T_t1:0), t ∈ { x, y, z }, for two The same position of width figure, remembers C_aNonzero eigenvalue number be M_a, remember C_bNonzero eigenvalue number be M_b, M_a、M_bAs tensor Rank of matrix it is approximate；

If 2.2.2, M_a=M_b, then two width figures have same type feature in this position, calculate the similarity of this position

${\mathrm{\γ}}_{ab}(x,y,z)=\frac{2\underset{p,q\∈\mathrm{\Ω}}{\Σ}\left|C_{a_H}\left(p\right)C_{b_H}^{*}\left(q\right)\right|}{\underset{p,q\∈\mathrm{\Ω}}{\Σ}\left|C_{a_H}\left(p\right)C_{a_H}^{*}\left(q\right)\right|+\underset{p,q\∈\mathrm{\Ω}}{\Σ}\left|C_{b_H}\left(p\right)C_{b_H}^{*}\left(q\right)\right|}---\left(3\right)$

Calculate threshold valueFusion rule is：

γ_abDuring≤α, this position is redundancy, selects weighted criterion：

C_fH=ω_aC_aH+ω_bC_bH (4)

${\mathrm{\ω}}_a=\frac{{\mathrm{\δ}}_a}{{\mathrm{\δ}}_a+{\mathrm{\δ}}_b},{\mathrm{\ω}}_b=\frac{{\mathrm{\δ}}_b}{{\mathrm{\δ}}_a+{\mathrm{\δ}}_b}$

${\mathrm{\δ}}_a={\mathrm{\σ}}_a\×\underset{x,y}{max}\left(abs\right(E_a\left(x,y,z\right)\left)\right),{\mathrm{\δ}}_b={\mathrm{\σ}}_b\×\underset{x,y}{max}\left(abs\right(E_b\left(x,y,z\right)\left)\right).$

γ_abDuring ＞ α, this position is complementary information, using MRE criterions:

$C_{fH}=\left\{\begin{array}{c}{C}_{aH},S_{aH}\&GreaterEqual;S_{bH}\\ C_{bH},S_{aH}<S_{bH}\end{array}\right.;---\left(5\right)$

$S_{tH}={\mathrm{\&sigma;}}_t=\frac{1}{N_{\mathrm{\&Omega;}}}\underset{p\&Element;\mathrm{\&Omega;}}{\&Sigma;}\sqrt{{\left(C_{tH}\right(p)-{\stackrel{\&OverBar;}{C}}_{tH})}^2},t\&Element;\{a,b\};$

If 2.2.3, M_a≠M_b, fusion criterion：

$C_{fH}=\left\{\begin{array}{c}{C}_{aH},M_a\&GreaterEqual;M_b\\ C_{bH},M_a<M_b\end{array}\right.---\left(6\right)$

5. it is according to claim 1 based on 3D shear wave conversion feature level Method of Medical Image Fusion, it is characterised in that The step 3 does backward shear conversion to the image after DWT or DTCWT inverse transformations, specific as follows：

For backward shear conversion refer to before to shear conversion inverse operation, wherein for z directions shear conversion refer to it is right Point in data carries out following coordinate transform：

It is for the shear transformation for mula in x directions：

Doing shear transformation for mula for y directions is：

Wherein, (x, y, z) is the coordinate before conversion, and (x ', y ', z ') is the coordinate after conversion.