PET and MR image fusion method

Abstract

The invention provides a PET (positron emission tomography) and MR (magnetic resonance) image fusion method, and belongs to the technical field of medical image processing. Carrying out image enhancement by using Radon transformation and a Gabor filter; performing signal decomposition on the enhanced image in a frequency domain through wavelet transform; establishing an image high-frequency component fusion model by adopting a pulse coupling neural network and calculating a high-frequency fusion sub-image of the image, and meanwhile, establishing an image low-frequency component fusion model by adopting a multi-dimensional region energy rule and calculating a low-frequency fusion sub-image of the image; performing wavelet inverse transformation on the fused image to complete image fusion; wherein the preprocessing comprises the steps of carrying out size normalization on the PET image and the MR image, carrying out bilinear interpolation to connect interlayer images, and carrying out gray feature normalization processing. According to the PET and MR image fusion method, the brightness of the fused image is improved, the image contrast is increased, the fusion details are improved, and the fusion effect is effectively improved.

Description

Technical Field

[0001] The invention belongs to the technical field of medical image processing, and in particular relates to a PET and MR image fusion method. Background Art

[0002] Image fusion is the process of combining two or more different images through algorithms to form a new image that meets certain conditions and describes the target or scene more accurately, comprehensively and reliably. Medical image fusion makes full use of the redundant and complementary information contained in the content, color, morphology and resolution of multiple original single-modality images to achieve complementary advantages. For example, the fusion of pseudo-color PET images and grayscale MR images can retain both the texture details of the soft tissue in the MR image and the functional color information provided by the PET image in the fused image, which is beneficial for the diagnosis of malignant diseases such as cancer and tumors.

[0003] For the original brain images collected, due to the differences in the scale, time, light, etc. of the captured medical images, the original images will have certain differences in time domain and space domain characteristics. Therefore, the reliability of the selected correlation features will ultimately affect the effect of image fusion.

[0004] For medical brain image research, texture features are very important. Therefore, image texture feature processing, image segmentation, and classification association fusion are very important. These image fusion methods can be roughly divided into the following categories: (1) image fusion algorithms based on spatial domain; (2) image fusion algorithms based on transform domain; (3) image fusion algorithms based on sparse representation; (4) image fusion algorithms based on deep learning, etc.

[0005] Wavelet transform is a commonly used method in image fusion. First, the wavelet transform algorithm is used to decompose and reconstruct the signal in the frequency domain of the source image. The features of the signal in each frequency band are extracted, and then the distribution energy of each frequency band obtained by calculation is used as the feature vector. Any image is decomposed into 4 parts after a single-layer two-dimensional discrete wavelet transform: 1 low-frequency part and 3 high-frequency parts. The low-frequency part mainly describes the outline of the image and concentrates most of the energy of the image. The high-frequency part is divided into components in three directions: high-frequency vertical direction, high-frequency horizontal direction, and high-frequency diagonal direction. The high-frequency part is mainly a detail signal. At the same time, wavelet decomposition is a multi-level decomposition, and the low-frequency part can be further decomposed into multiple layers. In image fusion based on wavelet transform, the processing method of the coefficients in the frequency band subband will greatly affect the effect of image fusion. Low frequency means that the color changes slowly, that is, the grayscale changes slowly, and it is a continuous gradient area. For the image of glioma, it is the brain, which is the background information. High frequency is an image with fast frequency changes, that is, the grayscale difference between adjacent areas is large. In the image of brain glioma, it is the edge and boundary between the lesion and the surrounding normal brain tissue. In addition, the lesion is different from the normal point in color, and the grayscale changes rapidly, which is the high-frequency part.

[0006] Pulse coupled neural network (PCNN), also known as the third generation artificial neural network, retains more detailed information. Its signal form and processing principle are more in line with the physiological characteristics of the human visual nervous system. It can be used in image processing fields such as image fusion, smoothing, segmentation, recognition, edge detection, and image fusion.

[0007] The present invention first normalizes the original image, then uses a median filter to perform filtering and noise reduction processing to retain more texture detail feature information, and then mainly decomposes the image into a series of sub-images in different scale spaces through wavelet transform, and decomposes the signal in the frequency domain of the source image. The characteristics of the signal in each frequency band are extracted, and then the distribution energy of each frequency band obtained by calculation is used as the feature vector to study its change characteristics at different scales. For low-frequency sub-images, the regional energy rule is adopted for fusion; for high-frequency sub-images, coefficient fusion is performed based on an adaptive pulse coupled neural network. Finally, the fused image is subjected to the inverse transform of the base wavelet transform to complete the image fusion. The fused medical image can be applied to the analysis, treatment and intraoperative navigation of diseases, and has a high research value. Summary of the invention

[0008] In order to realize the above-mentioned image fusion technology, the present invention firstly normalizes the original image, then uses a median filter to filter and reduce noise, retains more texture detail feature information, and then uses the improved wavelet transform to extract features from the normalized and filtered image group, mainly decomposing the image into a series of sub-images in different scale spaces, and decomposing and reconstructing the signal in the frequency domain of the source image. The features of the signal in each frequency band are extracted, and then the distribution energy of each frequency band obtained by calculation is used as the feature vector. Its change characteristics at different scales are studied. For low-frequency sub-images, the regional energy rule is used for fusion; for high-frequency sub-images, coefficient fusion is performed based on an adaptive pulse coupled neural network. Finally, the fused image is subjected to an inverse basic wavelet transform to complete the image fusion.

[0009] The present invention proposes a PET and MR image fusion method, comprising: performing pre-processing and median filtering on the PET and MR images; performing image enhancement by using Radon transform and Gabor filter; performing signal decomposition of the enhanced image in the frequency domain by using wavelet transform; establishing an image high-frequency component fusion model by using pulse coupled neural network, calculating the high-frequency fusion sub-image of the image, and at the same time, establishing an image low-frequency component fusion model by using multi-dimensional regional energy rule, calculating the low-frequency fusion sub-image of the image; performing base wavelet inverse transform on the fused image to complete image fusion; wherein the pre-processing comprises size normalization, bilinear interpolation connection of inter-layer images, and grayscale feature normalization processing on the PET and MR images.

[0010] Furthermore, the pulse coupled neural network reconstructs the high frequency components of the image through an adaptive parameter model, and the adaptive parameter model is:

[0011]

[0012] Among them, the set of neighborhood neurons X = |i,j|L i,j [n]>0∩|i,j|Y i,j [n] = 0, L i,j [n] is the output of the nth neuron coupling access field, Y ij [n] is the pulse sequence of the nth neuron coupling output, (i, j) is the coordinate of the pixel in the image group, I x is the neighboring neuron corresponding to the neuron in the input image I, S(i,j) is the stimulus factor, α is the control coefficient, synaptic connectivity coefficient β, m i (n) where i = 1 or 2, m1(n) and m2(n) represent the means corresponding to the ignition area and the unignition area, respectively, and the internal activity term U y [n] is the signal L y[n] is multiplied by the synaptic connectivity coefficient β, then added to the positive offset 1, and finally added to the signal F y [n] multiplied together.

[0013] Furthermore, in the image high-frequency component fusion model, after a neuron of the pulse-coupled neural network completes an output, the outputs of the neurons around the neuron are weighted summed to obtain a new connection domain, and the feedback domain and the connection domain are linearly multiplied to obtain an internal activity item; the internal activity item is compared with the dynamic threshold, and if the former is greater than the latter, the neuron generates a new output; this process is repeated until the pulse sequence no longer changes during the output; the threshold is updated, and a new round of iteration is started.

[0014] Furthermore, in the multi-dimensional regional energy rule to establish the image low-frequency component fusion model, the weighted fusion rule is:

[0015]

[0016] in, and are the k-th directional subband of the lth layer of the fused image F and the source images A and B respectively, and α(i,j) is the weighting coefficient.

[0017] The beneficial effects of this method are as follows: 1. The fusion rule of regional energy can better improve the brightness of the fused image than the traditional fusion rule of taking the average value of low frequency and taking the maximum value of high frequency, and the fusion rule of using different window functions can better obtain the information in the original image than the fusion rule of using the mean window function, improve the details, increase the image contrast, and effectively improve the fusion effect. 2. In order to better evaluate the fusion results, the objective evaluation methods of mean, variance, entropy and average cross entropy are used to improve the reliability of the system and the utilization rate of image information, reduce the complexity of the algorithm, reduce the time cost, and improve the effect and performance of the algorithm, which is real-time and universal. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] Figure 1 This is the overall flow chart of this patent.

[0019] Figure 2 These are the original PET and MR images of glioblastoma (the left image is the MR image, and the right image is the PET image).

[0020] Figure 3 These are images after median filtering (the left image is an MR image, and the right image is a PET image).

[0021] Figure 4 It is the Radon transform sinogram of the MR image and the enhanced image.

[0022] Figure 5It is the Radon transform sinogram of the PET image and the enhanced image.

[0023] Figure 6 is the wavelet decomposition of MR images.

[0024] Figure 7 It is the wavelet decomposition of PET image.

[0025] Figure 8 It is a comparison of 6 fusion algorithms (a: logical filtering method, b: weighted average method based on artificial fish school, c: image fusion based on twin convolutional neural network, d: image fusion based on deep stacked convolutional neural network, e: image fusion based on pulse coupled neural network, f: method of the present invention).

[0026] Fig. 9 It is an evaluation index of the experimental results of 6 fusion methods (Exp 1: logical filtering method, Exp 2: weighted average method based on artificial fish school, Exp 3: image fusion based on twin convolutional neural network, Exp 4: image fusion based on deep stacked convolutional neural network, Exp 5: image fusion based on pulse coupled neural network, Exp 6: the method of the present invention). DETAILED DESCRIPTION

[0027] The present invention is described in detail below based on the accompanying drawings and examples, but the specific implementation methods of the present invention are not limited thereto.

[0028] like Figure 1 The figure shows the overall flow chart of the method of the present invention, a PET and MR image fusion method, comprising:

[0029] Step 1: Image standardization preprocessing, according to the differences of the input original image, normalize the size and grayscale parameters. Specifically:

[0030] 1.1 Size normalization.

[0031] Due to differences in the scale, time, light, etc. of the medical images taken, the original images will have certain differences in time and space domain characteristics, and the proportion of the lesion area will often be different. Therefore, the reliability of the selected correlation features will ultimately affect the effect of image fusion. Then the overall image contrast effect will not be very good. Therefore, in order to facilitate subsequent processing, the size and grayscale of the brain lesion area images must be normalized during image preprocessing. The main purpose is to adjust the size and grayscale features of the collected brain images to a predetermined range, so that the input original images have a relatively close mean difference value, thereby ensuring that there is a same standard during subsequent processing.

[0032] Preferably, the size of the initial image is set to 128*128, and the size of the image that does not meet the standard will be scaled. This image size normalization process is actually a geometric transformation of the image.

[0033] 1.2 Use bilinear interpolation to connect images between layers.

[0034] After obtaining the original image with similar initial features, the image will be used as the original reference image, and the benchmark will also change. The target pixel calculated by this method can ensure that the entire target image has a relatively high feature correlation. However, the address of the pixel of the original image obtained may no longer be a positive number but a fraction, that is, the calculated original pixel coordinate corresponding to the target pixel may be a decimal, which is between the four adjacent pixels. Therefore, it is necessary to introduce the concept of grayscale interpolation to compensate for the error caused by image scaling.

[0035] In principle, the bilinear interpolation method is actually an improvement on the nearest neighbor method, which mainly interpolates the grayscale values ​​of the decimal coordinate point in four directions to calculate the value of the point. For the target image to be processed, first set a given feature parameter floating point coordinate to (i+m,j+n), where i and j are non-negative numbers, and m and n are coordinate floating point numbers. Then the pixel parameter L(i+m,j+n) here can be calculated by the four adjacent coordinates of the original image.

[0036] The coordinates of the four points are:

[0037]

[0038] The final value of the floating point coordinate will also be determined by the four adjacent pixel parameters:

[0039] L(i+m,j+n)=(1-m)(1-n)L(i,j)+(1-m)·n·L(i,j+1)+(1-n)·m·L(i+1,j)

[0040] +m·n·L(i+1,m+1)+G i,j (4)

[0041] Where L(i,j) represents the pixel value of the original image. i,j It represents the compensation factor, which is the evaluation value given by the fit of the given feature parameter at the interpolation point.

[0042] Although the computational complexity of the bilinear interpolation method is relatively high, in terms of medical image accuracy, this method has higher image quality when scaling, and the original feature parameters are used as an evaluation reference, so the normalized image obtained will not have pixel discontinuities.

[0043] In addition, the bilinear interpolation method is actually a filter. Due to the characteristics of low-pass filtering, it will form an over-filtering relative to the high-frequency signal part, and the image contour will appear blurred, so a certain compensation is given to the high-frequency part through the compensation factor. At the same time, the bilinear interpolation method is only used as an initial normalization of the image size for the entire system, and it does not have a very obvious impact on the grayscale characteristics of the brain lesion area, so this method has actually met the needs.

[0044] 1.3 Grayscale feature normalization

[0045] Since the wavelet transform requires the initial input image to be grayscale, the image also needs to be normalized after the image size is normalized. Due to the differences in the time when the image is taken, there are certain differences between different groups of images. If they are not processed, it will be very difficult and distorted to process different groups of images. The same principle as the image size normalization processing, the input original image needs to be converted into a normalized image with the same mean and variance.

[0046] The normalization function is:

[0047]

[0048] Among them, L(i+m,j+n) is the grayscale value before image processing, δ is the mean value before image processing, and μ is the grayscale variance before image processing. Then P(i+m,j+n) is the grayscale value of the image after normalization processing, δ0 is the grayscale mean value after processing, and μ0 is the grayscale variance value after processing.

[0049] After the size and grayscale value normalization of the above images, not only the size and grayscale mean are the same, but the overall image will also be clearer, which is conducive to subsequent fusion processing.

[0050] Step 2: Use a median filter to filter the normalized original image to retain more detail information.

[0051] Step 2: Median filter processing: For the filtering of the original image, considering that the final image needs to retain detail information, the median filter is used. The median filter is a nonlinear signal processing method. Under certain constraints, it can avoid the blurring of details caused by the linear filter, and it has better denoising effects for pulse interference and image scanning. The noise variance σ of the median filter output 2 ,for:

[0052]

[0053] where σ i 2 is the variance parameter of the input noise, d is the window size of the median filter, is the mean value of the input noise, is the density function of the input noise. It can be seen that the median filter output has a great relationship with the input noise density distribution function. The median filter is not very effective in dealing with random noise, but it has a good denoising effect on pulse interference, especially for narrow pulse interference within the window size range.

[0054] Step 3 uses Radon transform and Gabor filter to enhance the image.

[0055] Radon transform is a transformation to find the projection of a given image in the direction of a ray at a specified angle. For the projection of a two-dimensional image in the direction of a ray at a specified angle, it is actually the linear integral value in a certain direction. For a digital image, it is the cumulative sum in that direction.

[0056] There is a Radon transform of the function at any angle. Then the projection function of the Radon transform at any angle θ on the image can be expressed as:

[0057]

[0058] f(i ′ ,j ′ )=f(i ′ cosθ-j ′ sinθ,i ′ sinθ-j ′ cosθ) (8)

[0059] Where (i ′ ,j ′ ) is the result of rotating the image at (i, j) coordinates by an angle of θ, which means:

[0060]

[0061] Then in two-dimensional space, the Radon transformation function, that is, the projection at multiple angles in two-dimensional space, can be expressed as:

[0062] C θ,d =∫∫ D f(i,j)ξ(υ-icosθ-j sinθ)dυdθ (10)

[0063] Where D is the two-dimensional spatial domain representing the image, f(i,j) is the grayscale parameter of the image at position (i,j), ξ represents the impulse function. υ is the distance on the straight line, and after projection at an angle of θ, υ is i cosθ+j sinθ.

[0064] At the same time, for different pixel points (i, j) on a straight line, after Radon transformation, under given θ, υ conditions, they intersect at The sine curve, so as long as we determine You can get straight line information.

[0065] When an image undergoes geometric transformations in two-dimensional space, such as scaling, rotation, and translation, the Radon transform will maintain certain characteristics. For example, in the case of scaling:

[0066]

[0067] Where η is the variation parameter and α is the proportional factor

[0068] In the case of rotation, assume that the sequence of the original image after projection is C θ (k), then within a fixed period, the image f(i,j,b) after rotation by b is:

[0069] C[f(i,j,b)]=C θ+b (k) (12)

[0070] When inputting an image, we first need to perform a discrete Radon transform on the original image, then its n-th order moment is expressed as:

[0071] ψ n (θ)=∑ k k n υ θ (k) (13)

[0072] where ψ n (θ) is the nth order moment of the image, which is regarded as the quality of the image. θ (k) is the projection sequence of the image. θ (k) The centroid of the projection sequence can be expressed as:

[0073]

[0074] Then, the ordinary N-order moment after this basic transformation is not invariant to rotation, scaling, and translation. When the image is translated or rotated, the projection will also change accordingly. However, the central moment is invariant. In order to achieve invariance in the discrete Radon transform process, the ordinary moment needs to be replaced by the central moment to obtain invariance.

[0075]

[0076] in Yes θ (k) is the center of mass, then the central moment of order n is translationally invariant and is a constant.

[0077] Then when the image changes in proportion, the corresponding projection υ θ (k) will also change. Assume that the scaling factor is α and:

[0078]

[0079] but:

[0080] ψ n (θ) α =∑ k k n θ υ θ (k) α =α n+2 ψ n (θ) (17)

[0081] but:

[0082] ψ0(θ) α =α 2 ψ0(θ) (18)

[0083] ψ1(θ) α =α 3 ψ1(θ) (19)

[0084] Then ψ n (θ) α The center distance is:

[0085]

[0086] Then assume that the scale and translation invariant feature function

[0087] By substituting into the above formula, we can get:

[0088]

[0089] but get

[0090] So we get the scaling and translation invariant function

[0091]

[0092] Therefore, for Radon transform, the invariance result is achieved by constructing invariant moments.

[0093] The image after Radon transformation is filtered and enhanced using Gabor filter. Gabor filter has good frequency selectivity and direction selectivity, which is consistent with the characteristics of fingerprint image. Moreover, according to the characteristics of Gabor function that it is bandpass in the x direction and low pass in the y direction, rotating the Gabor filter can effectively remove the noise in the signal and achieve image enhancement.

[0094] The expression of the even-symmetric Gabor function is as follows:

[0095]

[0096] in:

[0097]

[0098] α is the direction of the Gabor filter, which is actually taken as the estimated image block direction; (x, y) is the coordinate of the pixel point; f is the image frequency, They are the spatial constants of the Gaussian envelope along the x and y axes, respectively, which are given by experience and are generally 4. The Radon transformed image is filtered by a Gabor filter, and then the filtered and enhanced CT and PET images are obtained by using the inverse Radon transform.

[0099] Step 4: Decompose the enhanced image into signals in the frequency domain through wavelet transform.

[0100] In step 4, the wavelet decomposition image is divided into 1 low-frequency sub-image and 3 high-frequency sub-images.

[0101] The explanation of wavelet transform on the function is that if the function f(t) can be expressed as a linear decomposition:

[0102] f(t)=∑ l a l ψ l (t) where: l is an integer index of a finite or infinite sum; a l is the real-valued expansion coefficient; ψ l (t) is the set of real-valued functions of parameter t, and is called ψ l (t) is the expansion set. If the above formula is unique, then this set can be called a basis of the function class that can be expanded. If the basis is orthogonal, that is:

[0103] < k (t),ψ l (t)>=∫ψ k (t)ψ l (t)dt=0(k≠l) (26)

[0104] Then, the coefficients can be calculated using the inner product:

[0105] a k = <f(t),ψ l (t)>=∫f(t)ψ k (t)dt (27)

[0106] For the Fourier series, the orthogonal basis functions ψ l (t) is sin(kω0t) and cos(kω0t) with kω0. For the wavelet expansion, a two-parameter system is constructed so that equation (27) becomes f(t) = ∑ k ∑ j a j,k ψ j,t (t), where j, k are integer indices, and ψ j,t (t) Wavelet expansion functions that usually form an orthogonal basis.

[0107] Step 5: Calculate the high-frequency fusion sub-image and the low-frequency fusion sub-image of the image; specifically, it includes the following two parallel steps.

[0108] Step 5.1: Use the pulse coupled neural network (PNCC) fusion model to obtain the high-frequency fusion sub-image.

[0109] Traditional image fusion methods mostly use pixel-level or local features as conditions, and then randomly select and judge fusion. The disadvantage is that it cannot take global features into account. In order to change this drawback, this patent reconstructs high-frequency images through a pulse-coupled neural network with adaptive connection coefficients. The pulse-coupled neural network tracks the fusion process, fully utilizes the characteristics of global coupling, performs adaptive intelligent analysis and judgment on the original image, and performs fusion processing.

[0110] Each neuron of PCNN is usually composed of three domains: coupling access domain, connection input domain and pulse generation domain. It is composed of two functional units: feedback input and connection input, and is connected to its adjacent neurons using synaptic connection weights. Each functional unit must perform iterative operations during operation and decreases exponentially during the iteration process.

[0111] For the feature data extracted in the early stage, the corresponding feature intensity and energy are calculated through the feature parameters. The feature intensity is expressed as:

[0112] P(i,j)=G(g(i,j)-g(i+a,j+b)) 2 +(g(i+a,j)-g(i,j+b)) (28)

[0113] Where (i, j) is the coordinate of the pixel in the image group, and g(i, j) is the grayscale parameter of the pixel at the corresponding coordinate (i, j). Then, while removing the corresponding feature intensity, the energy parameter of the feature can also be obtained. It can be expressed as:

[0114]

[0115] Where R(i,j) is the energy value, α i,j is the control factor. Since the differences in the fusion image group will cause certain imbalances after processing, in order to reduce the errors caused by the initial processing, a control factor is added when calculating the energy, and the adjustment is made according to the feature comparison similarity. Then the external stimulus input is:

[0116] S ij =P(i,j)+R(i,j) (30)

[0117]

[0118] Input F to the neuron feedback channel ij [n] can be expressed as:

[0119]

[0120] in, Represents the decay time constant of the input signal; V F represents the intrinsic potential of the neuron feedback channel; M ijkl Represents the access domain connection weight coefficient, V kl [n-1] represents the previous output of the neuron adjacent to the current neuron; S ij Represents the external input stimulus, which is the grayscale of the pixel f(i, j) corresponding to the pixel point of the input PET and MR images.

[0121] The output L of the neuron coupling access domain ij It can be expressed as:

[0122]

[0123] Among them, W ijkl Represents the decay time constant of the input signal; V L represents the intrinsic potential of the neuron feedback channel; in a single PCNN neuron, the signal F ij [n] and L ij[n] multiply and modulate, so that neurons can have strong connection weights between each other, and they do not necessarily affect each other. When there is a weak connection weight between them, neurons may affect each other. If one of two connected neurons does not work for a long time, even if there is a strong connection weight relationship between them, it is difficult for the other neuron to be stimulated and activated. If the connection weight between two neurons is relatively small, but one of them is constantly stimulated, the other neuron may also be activated.

[0124] Internal Activity ItemU ij [n] is the signal L ij [n] is multiplied by the synaptic connectivity coefficient β, then added to the positive offset 1, and finally added to the signal F ij [N] multiply. ij [N] is calculated as

[0125] U ij [N]=F ij [N]{1+βL ij [n]} (35)

[0126]

[0127] When the neuron connects the input domain to generate U ij When [n] is greater than the dynamic threshold λ, an output pulse sequence Y will be generated. ij [n].

[0128] The calculation formula of threshold λ is:

[0129]

[0130] Among them, V λ Usually we need to select the ratio U ij [n] An order of magnitude larger.

[0131] Y ij The calculation formula for [n] is:

[0132]

[0133] λ will gradually weaken according to the exponential law. When it weakens to U>λ, the neuron will be excited again, the dynamic threshold will increase again, and then the λ value will begin to gradually weaken again.

[0134] When a neuron inputs a pulse, the threshold of the neuron will be rapidly increased through feedback. When the threshold of the neuron is exceeded, the pulse generation will be turned off and the pulse will stop being emitted. Then the threshold index will begin to decrease. When the threshold is low, the pulse generator is turned on, the neuron is ignited, and is in an activated state, outputting a pulse or a series of pulse sets. Then it can be expressed as:

[0135]

[0136] Y=Step(U j -λ j ) (40)

[0137] In the PCNN model, the pulse transmission of the connection domain has a delay. In order to reduce the impact of the delay, self-adaptation is performed. After a neuron completes an output, the model performs a weighted summation of the outputs of the neurons around the neuron to obtain a new connection domain, and then linearly multiplies the feedback domain and the connection domain to obtain the internal activity term, and then compares the internal activity term with the dynamic threshold. If the former is greater than the latter, the neuron generates a new output. Repeat this process until Y no longer changes, then update the threshold and start a new round of iteration.

[0138] There are many parameters in the PCNN model. The current theory cannot clearly define the effect of changing each parameter. In addition, the commonly used parameters are also designed by human experience. Therefore, in the process of setting, it is inevitable to encounter a situation of constant debugging. At present, in the field of adaptive models, there are mainly two schools, one is parameter adaptation and the other is network iteration adaptation.

[0139] After the introduction of the model, the firing of neurons in the original simplified PCNN model needs to be completed in three steps, but now it can be completed in just one step. Neurons can stimulate surrounding neurons and send pulses in the same iteration process.

[0140] The PCNN model that adjusts the internal features of neurons through an adaptive parameter model is used to fuse the high-frequency components of the image. The parameter port of the traditional PCNN algorithm is generally a fixed empirical value, which limits the adaptability of PCNN in image fusion and also affects the fusion effect. Therefore, the present invention sets an adaptive parameter model:

[0141]

[0142] Where X is the set of neighborhood neurons, that is, X = |i,j|L i,j [n]>0∩|i,j|Y i,j [n] = 0. xis the neighboring neuron corresponding to the neuron in the input image I. S(i,j) is the stimulus factor, α is the control coefficient, m1(n) and m2(n) represent the means corresponding to the ignition area and the unignition area, respectively.

[0143] m c (n) = ∑ i,j∈Φ F i.j C=1,2 (42) where Φ1=|i,j|Y i.j [n]=0;Φ2=|i,j|Y i.j [n-1]=1.

[0144] Through the mean discrete optimization constraint, the optimal solution of β is obtained. The adaptive fusion strategy of PCNN can be realized by adjusting the internal feature U of the neuron. The PCNN with the above adaptive parameters is applied to the fusion of high-frequency components. The absolute value image of the high-frequency component is used as the input image.

[0145] The trigger time is updated each time through the time iteration formula:

[0146] T x,y [n] = T x,y [n-1]+Y x,y [n] (43)

[0147] At the same time, the image fusion subband is determined by the larger trigger time output. The present invention selects the larger value at the time of triggering. As high frequency component:

[0148]

[0149] In the formula, and Represent the high frequency subbands of image A and image B respectively, so the corresponding trigger times are and

[0150] Step 5.2: Using the multi-dimensional regional energy fusion model to fuse the low-frequency sub-images of the image;

[0151] The wavelet coefficients of the low-frequency sub-image are fused using the multi-dimensional region energy rule to obtain the low-frequency fused sub-image. First, the energy in the multi-dimensional region needs to be calculated separately, so the directional region energy can be expressed as:

[0152]

[0153] Where k represents the respective subband.

[0154] At the same time, the directional regional gradient and variance can be expressed as:

[0155]

[0156] Where μ is the average value of the coefficients in template A.

[0157] For the low-frequency subband, all coefficients in the subband are traversed, and the gradient values ​​of the eight neighborhoods of the two low-frequency subbands centered at (x, y) are calculated for each position (x, y). Then, the respective gradient values ​​are normalized, and the normalized gradient values ​​are used to perform coefficient weighting to determine the low-frequency fusion coefficient at (x, y);

[0158] Therefore, the low-frequency image fusion is performed by adopting the maximum energy selection principle. At the same time, each iteration uses the spatial frequency (SF) to improve the regional energy and The specific formula is

[0159]

[0160] Then the specific low-frequency coefficient fusion rules are as follows:

[0161]

[0162] Where X is the enhanced image, L X is the low-frequency coefficient matrix of the enhanced image after decomposition, and K is the number of wavelet decomposition layers.

[0163] We can obtain the values ​​at pixel (m,n) respectively. and The regional energy ratio of the low-frequency coefficient can be calculated as (m,n), as shown in the formula,

[0164]

[0165] The low-frequency coefficient matrix of image fusion is L F

[0166]

[0167] Where mean[·] is the calculated mean; std[·] is the calculated standard deviation.

[0168] The energy of the related area of ​​two images A and B is

[0169] E A,B (i,j)=∑ m ∑ n ω(m,n)I A (i+m,j+n)I B (i+m,j+n) (55)

[0170] Where: I A (i, j) and I B(i, j) are images A and B at coordinates I A (i, j) is the pixel value at position; ω(m, n) is the corresponding window function. When A and B are the same, E A,B It is the regional energy of a single image.

[0171] For the low-frequency part, a window function ω3 is used, in which the values ​​of the eight neighborhoods are all non-zero and the values ​​of the center and the neighborhood do not change much, to better match the WBCT neighborhood with relatively less drastic changes, because the low-frequency WBCT coefficients mainly contain lower-frequency information. ω1, ω2, and ω3 are

[0172]

[0173] Define the regional matching degree of two images as

[0174]

[0175] M A,B It reflects the approximation of the energy of the two image regions. When A and B are close, M A,B →1, and when A and B are quite different, M A,B →0. According to E A,B and M A,B , we get the weighted fusion rule

[0176]

[0177] Where: and are the k-th directional subbands of the lth layer of the fused image F and the enhanced images A and B respectively, α(i,j) is the weighting coefficient,

[0178]

[0179] Where T is A,B The relevant threshold is generally between 0.50 and 1.00.

[0180] Step 6: Perform inverse wavelet transform on the fused image to complete image fusion.

[0181] The high-frequency fusion sub-image and the low-frequency fusion sub-image obtained in step 5 are fused by inverse wavelet transform to construct the final fusion image.

[0182] The inverse transformation reconstruction of the two-dimensional Mallat algorithm image can be expressed as:

[0183]

[0184] In the formula, is a low-pass synthesis filter, It is a high-pass integrated filter.

[0185] So far, the present invention has completed the introduction of a PET and MR image fusion method. The following provides the experimental results of the method shown in the present invention, wherein Figure 2 The original PET and MR images of glioblastoma are shown (the left image is MR image, the right image is PET image). Figure 3 is the image after median filtering (the left image is MR image, the right image is PET image), Figure 4 is the Radon transform sinogram of the MR image and the enhanced image, Figure 5 is the Radon transform sinogram of the PET image and the enhanced image, Figure 6 is the wavelet decomposition of the MR image, Figure 7 is the wavelet decomposition of the PET image, Figure 8 The figure is a comparison between the method of the present invention and the other five common methods in the art (a: logical filtering method, b: weighted average method based on artificial fish school, c: image fusion based on twin convolutional neural network, d: image fusion based on deep stacked convolutional neural network, e: image fusion based on pulse coupled neural network, f: the method of the present invention). Fig. 9 The experimental results and evaluation indicators of these comparison algorithms are provided. It can be found that the experimental indicators of the improved algorithms of the present invention are all the best, which further proves the effectiveness of the improvements.

Claims

1. A PET and MR image fusion method, comprising: Performing median filtering on the PET and MR images after preprocessing; Image enhancement using Radon transform and Gabor filter; The enhanced image is decomposed in the frequency domain by wavelet transform; a pulse coupled neural network is used to establish a high-frequency component fusion model of the image, and the high-frequency fusion sub-image of the image is calculated. At the same time, a multi-dimensional regional energy rule is used to establish a low-frequency component fusion model of the image, and the low-frequency fusion sub-image of the image is calculated. The fused image is subjected to inverse wavelet transform to complete the image fusion; wherein the preprocessing includes size normalization of the PET and MR images, bilinear interpolation to connect interlayer images, and grayscale feature normalization.

2. The PET and MR image fusion method according to claim 1, characterized in that: The size normalization is to perform geometric transformation on the PET and MR images to obtain images with a size of 128*128.

3. The PET and MR image fusion method according to claim 1, characterized in that: The bilinear interpolation connection between the layer images is to calculate the pixel value of the image at the floating point coordinate pixel point (i+m, jj+n) using the following formula: L(i+m, jj+n)=(1-m)(1-n)L(i, j)+(1-m)·n·L(i, j+1)+(1-n)·m·L(i+1, j)+m·n·L(i+1, m+1)+G i，j , where L(i, j) represents the pixel value of the original image at coordinate (i, j), G i，j Represents the compensation factor, which is determined by the evaluation value of the fit of the given feature parameter at the interpolation point.

4. The PET and MR image fusion method according to claim 1, characterized in that: The grayscale feature normalization setting is: when the grayscale value L(i+m,γ+n) before image processing is greater than or equal to the grayscale mean value δ before image processing, the grayscale value of the image after grayscale normalization is When the gray value L(i+m, j+n) before image processing is less than the gray mean δ before image processing, the gray value of the image after gray normalization is Among them, μ is the gray value variance before image processing, and μ0 is the gray value variance after processing.

5. The PET and MR image fusion method according to claim 1, characterized in that: The noise variance σ output in the median filter is 2 for Among them, σ i 2 is the variance parameter of the input noise, and d is the filter window size of the median filter.

6. The PET and MR image fusion method according to claim 1, characterized in that: The Gabor filter has a Gaussian envelope with a spatial constant of 4 along the x and y axes.

7. The PET and MR image fusion method according to claim 1, characterized in that: The wavelet transform decomposes each input image into one low-frequency sub-image and three high-frequency sub-images.

8. The PET and MR image fusion method according to any one of claims 1 to 7, characterized in that: The pulse coupled neural network reconstructs the high frequency components of the image through an adaptive parameter model, and the adaptive parameter model is: Among them, the set of neighborhood neurons X = |i, j|L i，j [n]>0∩|i,j|Y i，j [n] = 0, L i，j [n] is the output of the nth neuron coupling access field, Y ij [n] is the pulse sequence of the nth neuron coupling output, (i, j) is the coordinate of the pixel in the image group, I x is the input image, the neighboring neurons corresponding to the neurons in the middle, S(i, j) is the stimulus factor, α is the control coefficient, β is the synaptic connectivity coefficient, m i (n) where i = 1 or 2, m1(n) and m2(n) represent the means corresponding to the ignition area and the unignition area, respectively, and the internal activity term U y [n] is the signal L y [n] is multiplied by the synaptic connectivity coefficient β, then added to the positive offset 1, and finally added to the signal F y [n] multiplied together.

9. The PET and MR image fusion method according to any one of claims 1 to 7, characterized in that: In the image high-frequency component fusion model, after a neuron of the pulse coupled neural network completes an output, the outputs of the neurons around the neuron are weighted summed to obtain a new connection domain, and the feedback domain and the connection domain are linearly multiplied to obtain an internal activity item; the internal activity item is compared with the dynamic threshold, and if the former is greater than the latter, the neuron generates a new output; this process is repeated until the pulse sequence no longer changes during output; Update the threshold and start a new round of iteration.

10. The PET and MR image fusion method according to any one of claims 1 to 7, characterized in that: In the multi-dimensional regional energy rule to establish the image low-frequency component fusion model, the weighted fusion rule is: in, and are the k-th directional subband of the lth layer of the fused image F and the source images A and B respectively, and α(i, j) is the weighting coefficient.

Images