A method for PET image harmonization based on image feature extraction

Abstract

The application provides a PET image and harmonization method based on image feature extraction, comprising the following steps: acquiring multiple groups of PET images of different sources, and pre-processing the PET images; performing gamma distribution parameter estimation on the pre-processed PET images to obtain multiple groups of image numerical distribution features; performing spatial autoregressive model coefficient estimation on the gamma residuals of each group after extracting the distribution features to obtain multiple groups of image spatial neighborhood relations; and performing image simulation based on the estimated gamma distribution features and spatial neighborhood relations of each group to obtain multiple groups of simulated image data at a unified numerical level. The application can not only help PET image multi-center research to improve the utilization rate of existing PET image data, but also help solve problems encountered during the equipment update transition period in clinical applications.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of image processing, and particularly relates to a PET image harmonization method based on image feature extraction. BACKGROUND

[0002] Positron emission tomography (PET) is a functional imaging method that can determine the metabolic level of tissues or organs by the enrichment level of a tracer in the human body. It has been widely used in the diagnosis, grading and prognosis of various types of tumor patients and the study of the human brain. Due to the low resolution of PET imaging, PET imaging is usually combined with other imaging techniques to compensate for this deficiency, such as PET / CT and PET / MR. In theory, PET images can exhibit the metabolic characteristics of various tissues and organs by using appropriate tracers. Therefore, researchers have developed tracers for different diseases or different targets by combining different isotopes with chemical substances to seek specific display of lesion sites in images. At the same time, in order to improve the reconstruction efficiency and quality of PET images, more reconstruction methods have been developed and applied to new generation scanning devices, such as ordered subset expectation maximization (OSEM) and point-spread function (PSF).

[0003] As a quantitative standard for PET images, the standard uptake value (SUV) has been proven to have certain differences between images obtained from different scanning devices, tracers and reconstruction methods. In addition, there are also influences from scanning settings and other aspects. Such differences can make the comparability between data from different sources lower, and severely limit the data source range of multi-center joint research of PET images. Although the comparability between data can be improved by strictly following uniform operation procedures, this will greatly reduce the utilization rate of existing PET image data.

[0004] The role of the harmonization algorithm for PET images is to process PET images from different sources to improve the comparability between images. The harmonization algorithm can improve the utilization rate of existing PET image data, reduce the burden of multi-center joint research on data acquisition, and avoid SUV value errors caused by operation differences in the scanning process as much as possible. Currently, internationally used and widely used harmonization algorithms include EQPET proposed by Siemens and GI-PET proposed by Canon.

[0005] The existing harmonization algorithms have the following disadvantages:

[0006] Use limitations of existing algorithms

[0007] EQPET and GI-PET are commercial harmonization algorithms, and there are certain limitations in the application range of the algorithms. EQPET is a built-in algorithm in the PET scanner of Siemens, and it can only be modified by setting at the beginning of scanning. In addition, its algorithm processing process is included in the image reconstruction process, and there is no corresponding display. Therefore, EQPET is currently only applicable to the data obtained by the scanner designed by Siemens. GI-PET can be applied to most PET scanners, but it is limited to use in Japan, and it has stopped selling at present.

[0008] Data application limitations of existing algorithms

[0009] EQPET and GI-PET are both harmonization algorithms based on three-dimensional Gaussian filtering, and the harmonization is based on the SUV reference range proposed by the European Nuclear Medicine Association and the Japanese Nuclear Medicine Association as the final standard. Since the filtering algorithm can only reduce the SUV value of the image in one direction, the above existing algorithms are not applicable to the data set below the reference range, so that this part of the data cannot be utilized. SUMMARY

[0010] To solve the above technical problems, the present application provides a PET image harmonization method based on image feature extraction, which extracts the numerical distribution characteristics and spatial neighborhood relationship of multiple groups of images through gamma distribution parameter estimation and spatial autoregressive model estimation, and then simulates the simulated image data at a unified numerical level according to the above results. This method not only helps PET image multicenter research to improve the utilization rate of existing PET image data, but also helps to solve the problems encountered in the transition period of device update in clinical application.

[0011] To achieve the above purpose, the present application provides a PET image harmonization method based on image feature extraction, comprising the following steps:

[0012] Obtain multiple groups of PET images of different sources and pre-process them;

[0013] Perform gamma distribution parameter estimation on each group of pre-processed PET images to obtain the estimated gamma distribution characteristics of each group;

[0014] Perform spatial autoregressive model coefficient estimation on the gamma residual of each group after extracting the gamma distribution characteristics to obtain the spatial neighborhood relationship;

[0015] Based on the estimated gamma distribution characteristics of each group and the spatial neighborhood relationship, image simulation is performed to obtain a plurality of groups of simulated image data at a unified numerical level.

[0016] Optionally, the pre-processing of the PET image comprises SUV correction, region of interest delineation, data normalization and filtering of the PET image.

[0017] Optionally, the method for SUV correction of the PET image is:

[0018] Numerical correction and attenuation correction are performed on the PET image.

[0019] The SUV-corrected PET image is normalized based on body weight.

[0020] Optionally, the calculation formula for SUV correction of the PET image is:

[0021] SUV bw = (image pixel value x correction slope + correction intercept) x SUV correction factor,

[0022] wherein,

[0023] Optionally, the method for region of interest delineation of the PET image is:

[0024] The normalized PET image is converted into an elliptical cylindrical data.

[0025] Optionally, the method for data normalization of the PET image is:

[0026] The elliptical cylindrical data is converted into a form with an overall mean of one.

[0027] Optionally, the method for filtering of the PET image is:

[0028] The elliptical cylindrical data after data normalization is denoised using a two-dimensional gamma filtering algorithm.

[0029] Optionally, the maximum likelihood objective equation for gamma distribution parameter estimation of each group of the pre-processed PET image is:

[0030]

[0031] wherein, x ijk is the length n of the randomly generated gamma distribution data subject to each voxel, μ ijk and φ ijk are calculated according to the distribution subject to each voxel.

[0032] Optionally, the spatial autoregressive model coefficient estimation is performed on the gamma residual after the gamma distribution feature extraction of each group, and the method for obtaining the spatial neighborhood relationship comprises the following steps:

[0033] Based on the spatial autoregressive model, the linear relationship between the neighborhoods is obtained.

[0034] Based on the linear relationship between the neighborhoods, a likelihood target equation is established to estimate the coefficients of the corresponding neighborhoods.

[0035] The target equation is minimized by using a nonlinear weighted least squares estimation method to obtain a neighborhood relationship model containing the coefficients of the corresponding neighborhoods.

[0036] The neighborhood relationship model is applied to the gamma estimation residual to extract the corresponding neighborhood information and calculate the residual.

[0037] If the residual still shows residual information, the neighborhood relationship model is adjusted for re-estimation, and if the preset requirements are met, the final neighborhood relationship model is saved.

[0038] Based on the neighborhood relationship model meeting the preset requirements, the spatial neighborhood relationship is obtained.

[0039] Compared with the prior art, the present application has the following advantages and technical effects:

[0040] The present application is based on the image feature extraction and harmonization algorithm, so it does not need external reference range as the harmonization standard, and it also does not have the one-way harmonization limitation of EQPET and GI-PET. In PET multicenter research, the application of the present algorithm can harmonize PET data from multiple sources and convert them to a unified and comparable level for disease research, instead of being limited to data from a unified source and strictly following a unified operation process, thereby improving the utilization rate of existing PET data in multicenter research.

[0041] In current clinical applications, due to the update of in-hospital equipment, the SUV value of the lesion obtained by scanning the existing patients after treatment on the new PET scanning equipment is higher than the SUV value of the lesion obtained by scanning the existing patients before treatment on the old PET scanning equipment. Therefore, it is difficult for doctors to judge the treatment effect of patients according to the PET data of two scans, and it is difficult for doctors to make prognosis for patients. The present application converts the PET images of two scans to a unified numerical level through harmonization, reduces the SUV value difference caused by different factors such as scanning equipment, tracer and reconstruction method, makes the results of two scans comparable, and solves the problems encountered in the transition period of equipment update in clinical applications. BRIEF DESCRIPTION OF DRAWINGS

[0042] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings:

[0043] Figure 1 This is a schematic diagram of a PET image harmonization method based on image feature extraction according to Embodiment 1 of the present invention. Detailed Implementation

[0044] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0045] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0046] Example 1

[0047] like Figure 1 As shown, this invention provides a PET image harmonization method based on image feature extraction, comprising the following steps:

[0048] Multiple sets of PET images from different sources were obtained and preprocessed.

[0049] The gamma distribution parameters of each group of preprocessed PET images are estimated to obtain the estimated gamma distribution characteristics of each group.

[0050] Spatial autoregressive model coefficients are estimated for the gamma residuals after extracting gamma distribution characteristics of each group to obtain spatial neighborhood relationships.

[0051] Image simulations are performed based on the estimated gamma distribution characteristics of each group and the spatial neighborhood relationships to obtain multiple sets of simulated image data at a uniform numerical level.

[0052] This embodiment uses two PET scanners with different image numerical levels as an example and describes them separately in four main steps according to the processing order. It should be noted that this invention can harmonize data obtained from various different scanners; here, only two scanning devices are used as an example. Since the first three steps (preprocessing, gamma parameter estimation, and spatial autoregressive model coefficient estimation) are performed separately for each set of PET data, the descriptions of the steps will only mention one set of data.

[0053] Preprocessing is mainly divided into four parts: SUV correction, delineation of region of interest (ROI), data normalization, and filtering.

[0054] The SUV value standard used in this embodiment is the standard uptake value body weight (SUVbw). The purpose of SUV correction is to correct the image numerically and decay, and to standardize the value on the basis of body weight to ensure the comparability between data. The correction method of SUV value is calculated according to the following equation:

[0055] SUVbw= (image pixel value x correction slope + correction intercept) x SUV correction factor

[0056] SUVbw= (image pixel value x correction slope + correction intercept) x SUV correction factor

[0057] Drawing the ROI refers to drawing the human body region around the lesion with the same ellipse, converting the human body shape data in the PET image into the data of an elliptical cylinder and inputting into the subsequent calculation step. The purpose of this step is to exclude the influence of the region outside the human body in the image on the subsequent estimation, and to reduce the calculation amount of the algorithm by reducing the number of voxels to be estimated, thereby improving the efficiency of the algorithm.

[0058] The normalization processing of data is to convert the data in the elliptical cylinder ROI into a form with a total mean value of one to improve the convergence speed of the algorithm.

[0059] The filtering processing is to use a two-dimensional gamma filtering algorithm to reduce the noise of the data to reduce the influence of noise on the algorithm.

[0060] Gamma distribution parameter estimation

[0061] In this embodiment, the parameterized gamma distribution is considered to represent the reconstructed activity value of the PET image and it is assumed that each voxel is subject to a respective independent gamma distribution. Let {z ijk ,i=1,…,N; j=1,…,N; k=1,…,K} be the SUV value corresponding to each voxel. Where N is the length of the pixel of each image, and K is the number of images in the group. Each voxel is subject to an independent gamma distribution, that is:

[0062] E(z ijk ) = μ ijk , Var(z ijk ) = μ ijk φ ijk

[0063] Where μ ijk and φ ijkis calculated according to the distribution that each voxel is subjected to. For the parameter estimation of the gamma distribution of each point, the maximum gamma likelihood estimation method and Newton update process are selected to optimize the best parameters. The target optimization equation is obtained from the probability density function of the gamma distribution:

[0064]

[0065] where x ijk is the length of the randomly generated gamma distribution data of n that each voxel is subjected to.

[0066] The specific optimization process is as follows:

[0067] Randomly generate the gamma distribution that each voxel is subjected to

[0068] Assuming that {block} represents the target voxel and the surrounding adjacent circle of voxels, the mean {Mean(block)} and variance {Variance(block)} of {block} are calculated. The shape parameter and scale parameter of the gamma distribution that the target voxel is subjected to can be determined as {[Mean(block) 2 / Variance(block)} and {Variance(block) / Mean(block)} respectively. It should be noted that the gamma distribution calculated by {block} is theoretically the distribution that the mean of {block} is subjected to. In the subsequent steps, the SUV value simulated by the estimated information will be biased compared to the SUV value corrected by the image. Therefore, after generating the random distribution, the ratio of the mean value of {block} to the original value will be calculated and saved for subsequent data simulation.

[0069] Initialize the target parameters μ ijk and φ ijk

[0070] The initialization value of the target parameters is calculated according to the randomly generated gamma distribution that each voxel is subjected to {x ijk}. The initialization value of μ ijk is Mean(x ijk ), and the initialization value of φ ijk is Variance(x ijk ) / Mean(x ijk ). Since we have normalized the data, the initialization value of μ ijk should be around 1.

[0071] Minimize the target equation

[0072] The implementation of the minimization objective function is based on the Newton update step, which iterates continuously by first updating μ and then updating φ until the parameter changes are below the tolerance. The first order {g(·)} and second order {h(·)} partial derivatives of the objective function with respect to μ and φ are shown as follows:

[0073]

[0074]

[0075]

[0076]

[0077] Taking the optimization of μ as an example, the Newton update step is calculated by where μ 1 is the updated parameter, μ 0 is the previous parameter, and λ is the step factor introduced between (0, 1]. The selection criterion of λ is to satisfy the minimum value of μ 1 greater than zero. The calculation of the tolerance is subject to When and only when the tolerances of μ and φ are both less than 10 -5 , the iteration stops.

[0078] In addition, the embodiment additionally adds a probability conversion step when calculating the residual error of the gamma distribution estimate to avoid the possible residual skewness. The calculation of the residual error is carried out according to the following equation:

[0079]

[0080] where Φ -1 (·) represents the inverse cumulative distribution function of the standard Gaussian distribution, and F(·) represents the cumulative distribution function of the gamma distribution to which each voxel is subject.

[0081] Spatial auto regressive (SAR) model parameter estimation

[0082] The embodiment analyzes the three-dimensional covariance of the gamma residual error through the spatial auto regressive model to obtain the residual spatial neighborhood relationship between the voxels and their neighborhoods. The SAR model specifies a linear relationship between a group of neighborhoods:

[0083]

[0084] where u(n) represents the gamma estimate residual error value of the target voxel, u(n-k) represents the adjacent voxel of the target voxel, and k=(k1, k2, k3). ∈(n) represents a group of variances σ2 of Gaussian white noise, while θ k denotes the coefficients of the corresponding neighborhood. For the neighborhood coefficients θ k of the target voxel, we consider a likelihood-based approach. The corresponding likelihood objective function is given by:

[0085]

[0086] where, denotes the voxels contained in the whole including the target voxel, is the voxels contained in the neighborhood excluding the target voxel. is a set of linear difference operators satisfying denotes the three-dimensional spectral density of the SAR process (in the following in the form of f θ (λ)), where P θ (λ) =∑ k θ k e iλ·k is the three-dimensional discrete Fourier transform of the SAR model coefficients. Taking the derivative of the objective function with respect to θ and setting the derivative to zero, we obtain:

[0087]

[0088] where, indexes the voxels of the non-zero coefficients. By Parseval’s relation, the equation can be transformed into the following form:

[0089]

[0090] where, is the sample estimate of the covariance, (l'-l|θ) =∫e -iλ·(l-l′) f θ (λ) dλ is the inverse Fourier transform of the spectral density f θ (λ), which gives the auto-covariance of the three-dimensional model.

[0091] The objective function is minimized by using a non-linear weighted least squares estimation to obtain the neighborhood relationship model containing the corresponding neighborhood coefficients. The obtained neighborhood model is applied to the gamma estimation residuals to extract their neighborhood relationships and calculate the residuals. The effectiveness of the neighborhood model is determined by residual analysis. If the residuals still show residual information, the neighborhood model is adjusted and re-estimated. If the requirements are met, it is saved as the final neighborhood model for data simulation.

[0092] Data simulation and numerical harmonization

[0093] The starting point of data simulation is Gaussian white noise, so the algorithm can realize the interval estimation of simulated image data by generating multiple sets of white noise for simulation at the same time. First, the neighborhood relationship is applied to the white noise:

[0094]

[0095] where, and respectively represent the fast Fourier transform (FFT) and its inverse transform (IFFT). The simulated data added with the neighborhood relationship is taken as the gamma estimation residual, and the estimated gamma distribution feature information is combined to simulate the normalized SUV value.

[0096] z ijk =F -1 (Φ(u(n))|μ ijk ,φ ijk )

[0097] where, F -1 (·) represents the inverse cumulative distribution function of the gamma distribution to which the voxels are subject. It should be noted that before applying the estimated gamma feature information, the mu ijk and phi ijk estimated from two sets of data of different sources need to be adjusted to approximate distribution by histogram matching method.

[0098] Restoring the simulated data from normalization requires additional multiplication by two coefficients. The first coefficient is the inverse of the proportion of the original value lost in preprocessing and distribution calculation, used to restore the numerical loss generated in the preprocessing process. The second coefficient is the mean of the divisors of the two data in normalization, used to restore it to a unified numerical level.

[0099] The two sets of data finally obtained will be at a unified level in numerical value, and at this time the numerical difference of the lesion area is considered to be the real difference between the lesions, significantly reducing the difference brought by the difference of scanning equipment, tracer, reconstruction method and other factors.

[0100] The present application provides a PET image harmonization method based on image feature extraction. The core idea of the present application is not to take external SUV reference range as the standard, but to harmonize the PET images obtained from different PET devices, tracers, reconstruction methods, etc. by an algorithm based on image feature extraction, so as to improve the comparability between data of different sources.

[0101] The method is based on an image feature information extraction algorithm, and a harmonized image at a unified numerical level is simulated through acquired feature information.

[0102] The above is only a preferred specific embodiment of the application, but the protection scope of the application is not limited thereto, and any person skilled in the art can easily think of changes or replacements within the technical range disclosed by the application, which should be covered in the protection scope of the application. Therefore, the protection scope of the application should be subject to the protection scope of the claims.

Claims

1. A method of PET image harmonization based on image feature extraction, characterized in that, Includes the following steps: Multiple sets of PET images from different sources were obtained and preprocessed. The gamma distribution parameters of each group of preprocessed PET images are estimated to obtain the estimated gamma distribution characteristics of each group. Spatial autoregressive model coefficients are estimated for the gamma residuals after extracting gamma distribution characteristics of each group to obtain spatial neighborhood relationships. Image simulation is performed based on the estimated gamma distribution characteristics of each group and the spatial neighborhood relationship to obtain multiple sets of simulated image data at a uniform numerical level. Consider using a parametric gamma distribution to represent the reconstructed activity values of PET images and assume that each voxel follows a respective independent gamma distribution; let {z ijk , i = 1,..., N; j = 1,..., N; k = 1,..., K} be the SUV value of each voxel; where N is the length of pixels of each image, K is the number of images in the group; each voxel follows an independent gamma distribution, i.e.: where μ ijk and φ ijk The parameters of the gamma distribution are calculated according to the distribution to which each voxel is subjected; for the estimation of the parameters of the gamma distribution for each point, the maximum gamma likelihood estimation method and the Newton update process are chosen to optimize the best parameters; the target optimization equation is obtained from the probability density function of the gamma distribution: where x ijk is the length of the randomly generated gamma distribution of length n that each voxel is subject to; The specific optimization process: Randomly generate the gamma distribution that each voxel follows. Assuming {block} represents the target voxel and its surrounding neighboring voxels, the mean {Mean(block)} and variance {Variance(block)} of {block} are calculated; the shape parameter and scale parameter of the gamma distribution that the target voxel obeys are determined as {[Mean(block)] 2 / Variance(block)} and {Variance(block) / Mean(block)} respectively; it is noted that the gamma distribution calculated by {block} is theoretically the distribution that the mean of {block} obeys; in the subsequent steps, the SUV value simulated by the estimated information will be deviated compared to the corrected SUV value of the image; after the random distribution is generated, the ratio of the mean and the original value of {block} will be calculated and saved for subsequent data simulation; Initialization of the target parameter μ ijk and φ ijk The initial value of the target parameter is calculated according to the random generated gamma distribution x ijk} to which each voxel is subjected; the initial value of μ ijk is Mean(x ijk ), and the initial value of φ ijk is Variance(x ijk ) / Mean(x ijk ); the data is normalized, and the initial value of μ ijk is around 1. Minimize the objective equation The minimization of the objective equation is achieved using an optimization algorithm based on the Newton update step. This involves iterating continuously in the order of first updating μ while keeping φ constant, then updating φ while keeping μ constant, until the parameter variation falls below the tolerance, at which point the optimal parameter selection is determined. The first-order {g(·)} and second-order {h(·)} partial derivatives of the objective equation with respect to μ and φ are shown below: For the case of μ-optimization, the Newton update step is computed by where μ 1 is the updated parameter, μ 0 is the previous parameter, and λ is a step factor introduced between (0, 1]. The criterion for choosing λ is to satisfy the minimum value of μ 1 greater than zero; the tolerance is computed according to The iteration stops when and only when the tolerances of μ and φ are both less than 10 -5 . An additional probability transformation step is added when calculating the residuals of the gamma distribution estimate. The residual calculation is as follows: where Φ -1 (·) denotes the inverse cumulative distribution function of the standard Gaussian distribution, F(·) denotes the cumulative distribution function of the Gamma distribution that each voxel is assumed to follow; Spatial Autoregressive Model Parameter Estimation The three-dimensional covariance of the gamma residuals is analyzed using a spatial autoregressive model to obtain the residual spatial neighborhood relationships between voxels and their neighbors; the SAR model specifies a set of linear relationships between neighborhoods: where u(n) represents the gamma estimation residual value of the target voxel, u(n-k) represents the neighboring voxel of the target voxel, and k = (k1, k2, k3); ∈(n) represents a set of Gaussian white noise with variance σ 2 , and θ k represents the coefficient of the corresponding neighborhood; for the estimation of the neighborhood coefficient θ k , a likelihood-based method is adopted; the corresponding likelihood objective equation of the method is: where N0 denotes the voxels contained in the whole including the target voxel, and N / N0 denotes the voxels contained in the neighborhood excluding the target voxel; P θ is a set of linear difference operators, satisfying P θ u(n) = ∈(n); is the three-dimensional spectral density of the SAR process, expressed in f θ is simplified in the form P θ (λ) = ∑ k θ k e iλ·k is the three-dimensional discrete Fourier transform of the SAR model coefficients; and the target equation is derived by taking the derivative of the equation with respect to θ and setting the derivative equal to zero, i.e. wherein, The corresponding voxel of the index nonzero coefficient; by Parseval's theorem, the equation is converted to: where is a sample estimate of the covariance, (l'-l|θ) = ∫e -iλ·(l-l') f θ (λ) dλ is the inverse Fourier transform of the spectral density f θ (λ) gives the autocovariance of the three-dimensional model; The objective equation is minimized using nonlinear weighted least squares estimation to obtain a neighborhood relationship model containing the corresponding neighborhood coefficients. The obtained neighborhood model is then applied to the gamma estimation residuals to extract the neighborhood relationships and calculate the residuals. The effectiveness of the neighborhood model is determined through residual analysis. If the residuals show that there is still residual information, the neighborhood model is adjusted and re-estimated. If the requirements are met, the model is saved as the final neighborhood model for data simulation. Data simulation and numerical harmonization The starting point for data simulation is Gaussian white noise. Multiple sets of white noise are generated simultaneously for simulation to achieve interval estimation of the simulated image data. First, neighborhood relationships are applied to the white noise: wherein, and respectively represent the fast Fourier transform FFT and the inverse transform IFFT; the analog data added with the neighborhood relationship is taken as the gamma estimation residual, and the normalized SUV value can be simulated by combining the estimated gamma distribution characteristic information. z ijk = F -1 (Φ(u(n))|μ ijk ,φ ijk ) where F -1 (·) denotes the inverse cumulative distribution function of the gamma distribution to which each voxel is subject to; it is to be noted that the μ ijk and φ ijk need to be adjusted to an approximate distribution using a histogram matching approach; To recover simulated data from normalization, two additional coefficients are needed. The first coefficient is the reciprocal of the proportion of loss in the original value during preprocessing and distribution calculation, used to recover the numerical loss incurred during preprocessing. The second coefficient is the mean of the divisors of the two data points during normalization, used to restore them to a uniform numerical level. The two sets of data obtained will be at the same numerical level. The numerical difference in the lesion area is considered to be the real difference between lesions, significantly reducing the differences caused by factors such as different scanning equipment, tracers, and reconstruction methods.

2. The method of claim 1, wherein, The preprocessing of the PET image includes: SUV correction, region of interest delineation, data normalization, and filtering.

3. The method of image feature extraction based PET image harmonization according to claim 2, characterized in that, The method for performing SUV correction on the PET image is as follows: Numerical up-correction and attenuation correction are performed on the PET image; The PET images were numerically and attenuated based on body weight.

4. The method of image feature extraction based PET image harmonization according to claim 3, characterized in that, The formula for calculating SUV correction of the PET image is as follows: SUVbw = (Image pixel value × Correction slope + Correction intercept) × SUV correction factor wherein 5. The image feature extraction based PET image harmonization method of claim 3, wherein, The method for delineating a region of interest on the PET image is: The normalized PET image is converted into an elliptical cylindrical data.

6. The image feature extraction based PET image harmonization method of claim 5, wherein, The method for data normalization processing on the PET image is: The elliptical cylindrical data is converted into a form with an overall mean of one.

7. The image feature extraction based PET image harmonization method of claim 6, wherein, The method for filtering processing on the PET image is: The elliptical cylindrical data after data normalization processing is denoised by using a two-dimensional gamma filtering algorithm.

8. The image feature extraction based PET image harmonization method of claim 1, wherein, The method for obtaining the spatial neighborhood relationship based on the spatial autoregressive model coefficient estimation of the gamma residual after extracting the gamma distribution characteristics of each group includes: Based on the spatial autoregressive model, a linear relationship between neighborhoods is obtained. Based on the linear relationship between the neighborhoods, a likelihood objective equation is established to estimate the coefficients of the corresponding neighborhoods. By using a nonlinear weighted least squares estimation method, the objective equation is minimized to obtain a neighborhood relationship model containing the coefficients of the corresponding neighborhoods. The neighborhood relationship model is applied to the gamma estimation residual to extract the corresponding neighborhood information and calculate the residual. If the residual still has residual information, the neighborhood relationship model is adjusted for re-estimation; if the preset requirements are met, the final neighborhood relationship model is saved. Based on the neighborhood relationship model meeting the preset requirements, the spatial neighborhood relationship is obtained.

Images