Magnetic particle imaging method based on multi-harmonic energy-full width at half maximum combined weighting

Abstract

The invention discloses a magnetic particle imaging method based on multi-harmonic energy-full width at half maximum combined weighting, which comprises the following steps: acquiring point spread functions of a plurality of harmonic waves with different frequencies of magnetic nanoparticles, and constructing a corresponding system matrix; measuring a reference sample which does not contain the magnetic nanoparticles, and obtaining a plurality of noise matrixes corresponding to harmonic waves with different frequencies; measuring the sample to be measured, and obtaining magnetization response voltage vectors of a plurality of harmonic waves with different frequencies; respectively calculating an energy two-norm of a system matrix, an energy two-norm of a noise matrix and a full width at half maximum of a point spread function of the system matrix and the noise matrix; splicing the magnetization response voltage signals of each harmonic system matrix to form a total system matrix and a total measurement signal vector; calculating a weight factor of each harmonic wave and forming a weight matrix; and constructing a solution equation of the particle concentration distribution of the sample to be detected, and carrying out iterative solution on the solution equation to obtain a particle concentration distribution image of the sample to be detected. According to the invention, efficient utilization of multi-harmonic information is realized.

Description

Technical Field

[0001] This invention relates to the field of magnetic nanoparticle imaging technology, and in particular to a magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting. Background Technology

[0002] Magnetic nanoparticle imaging (MPI) is a novel imaging technique based on the magnetization dynamics of superparamagnetic nanoparticles. This method applies an alternating magnetic field to induce significant nonlinear magnetization responses in superparamagnetic iron oxide nanoparticles (SPIONs) within localized unsaturated regions, and then reconstructs their spatial distribution using the voltage signal induced by a receiving coil. Since biological tissues themselves do not generate background signals, MPI enables direct, quantitative imaging of magnetic tracers. Compared to existing medical imaging techniques such as CT, MRI, PET, and ultrasound, MPI offers high detection sensitivity and real-time dynamic imaging capabilities, while avoiding limitations imposed by ionizing radiation and tissue attenuation effects. Therefore, MPI demonstrates significant technological advantages and application potential in areas such as non-invasive blood flow monitoring, micro-dose drug tracking, and cell and drug delivery carrier monitoring.

[0003] Currently, MPI imaging mainly includes two categories: x-space methods based on time-domain mapping and system matrix methods based on the frequency domain. The x-space method is based on the direct mapping between the time-domain signal and the FFP (field-free point) location, offering fast reconstruction speed, but it has high requirements for system linearity, scan trajectory, and magnetic field homogeneity. The system matrix method, by acquiring multiple harmonic components and constructing a pre-calibrated system matrix, achieves high-precision solutions for nanoparticle concentration distribution. Compared to the x-space method, the system matrix method offers higher spatial resolution and quantitative consistency, but requires higher calibration costs and introduces greater computational complexity.

[0004] The system matrix method is essentially a highly ill-conditioned linear inverse problem. Its singular values ​​decay rapidly, resulting in a large condition number for the system matrix, making it extremely sensitive to noise. Direct matrix inversion significantly amplifies higher-order harmonic noise, leading to reconstruction distortions such as fringe artifacts and blurred edges. Therefore, stabilization through regularization and iterative algorithms is necessary. Furthermore, the strong correlation between rows of the system matrix slows down the convergence speed of the iterative solver; simultaneously, the dominance of low-frequency components in the system matrix makes it difficult to effectively reconstruct high-frequency structural details, further limiting the improvement of imaging spatial resolution.

[0005] To overcome the aforementioned problems, existing research has begun to introduce weighting strategies into the system matrix reconstruction framework. By applying appropriate weights to different frequency components, the ill-conditioned nature caused by the rapid decay of singular values ​​in the system matrix can be mitigated to some extent, improving the stability of the inverse problem and effectively suppressing noise amplification effects. Furthermore, weighting can improve the sparsity and orthogonality of the system matrix, thereby accelerating the convergence speed of iterative reconstruction algorithms. Moreover, weighting can redistribute the relative contributions of each harmonic during the reconstruction process, appropriately increasing the weight of higher harmonics in the system matrix to enhance spatial resolution.

[0006] However, most existing weighting strategies only normalize the amplitude or energy of harmonics, failing to incorporate key imaging factors such as the spatial resolution characteristics and noise sensitivity of different harmonics into the weighting system. Since low-order harmonics have high energy but limited resolution, while high-order harmonics have high resolution but low signal-to-noise ratio, traditional energy-based weighting methods struggle to establish an effective dynamic balance between the two. Weighting based on the energy L2 norm leads to insufficient contribution from high-order harmonics, resulting in lower reconstructed image resolution; weighting based on the square of the energy L2 norm leads to excessive weighting of high-order harmonics, introducing numerous artifacts. Therefore, existing weighting methods still have significant limitations in improving the stability and imaging quality of the system matrix method. Summary of the Invention

[0007] To address the technical challenges in magnetic nanoparticle imaging, such as the difficulty in achieving both high spatial resolution and high signal-to-noise ratio, the strong ill-conditioning of the system matrix method, the easy amplification of noise, and the difficulty in effectively utilizing higher-order harmonics, this invention proposes a magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting. This method enables efficient utilization of multi-harmonic information, significantly improves the spatial resolution of reconstructed images, and suppresses artifacts.

[0008] To achieve the above objectives, this invention provides a magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting, comprising:

[0009] S1. In the magnetic particle imaging system, a magnetic nanoparticle sample placed at a zero magnetic field point is excited by an excitation magnetic field, and the zero magnetic field point is moved to traverse the imaging field of view. The point spread function of several harmonics of different frequencies of the magnetic nanoparticles is measured and obtained. Based on the point spread function of each harmonic, the corresponding system matrix is ​​constructed.

[0010] S2. Under the same imaging conditions as in S1, measure a reference sample without magnetic nanoparticles and obtain noise matrices corresponding to several different frequency harmonics;

[0011] S3. Under the same imaging conditions as S1, measure the sample to be tested and obtain the magnetization response voltage vectors of several different frequency harmonics.

[0012] S4. For each harmonic, calculate the energy norm of the system matrix, the energy norm of the noise matrix, and the full width at half maximum (FWHM) of its point spread function.

[0013] S5. The harmonic system matrices obtained in S1 are spliced ​​together with the magnetization response voltage signals obtained in S3 to form the total system matrix and the total measurement signal vector; and based on the energy norm 2 of the system matrix of each harmonic, the energy norm 2 of the noise matrix, and the full width at half maximum (FWHM) of the system matrix, the weighting factor of each harmonic is calculated and a weighting matrix is ​​constructed.

[0014] S6. Using the weight matrix, the total system matrix, and the total measurement signal vector, construct the solution equation for the particle concentration distribution of the sample to be tested, and iteratively solve the solution equation to obtain the particle concentration distribution image of the sample to be tested.

[0015] Preferably, the point spread function of the magnetic nanoparticles for several different frequency harmonics includes the point spread function of the third harmonic, fifth harmonic, seventh harmonic and ninth harmonic; the noise matrix corresponding to the several different frequency harmonics includes the noise matrix of the third harmonic, fifth harmonic, seventh harmonic and ninth harmonic.

[0016] Preferably, the overall system matrix is ​​constructed by concatenating the system matrices of each harmonic row by row, specifically as follows:

[0017] ;

[0018] In the formula, The system matrix for the third harmonic of magnetic nanoparticles. The system matrix is ​​for the fifth harmonic. The system matrix for the seventh harmonic is... The system matrix is ​​for the ninth harmonic. This is the overall system matrix;

[0019] The total measurement signal vector is constructed by concatenating the magnetization response voltage signal vectors of each harmonic row by row, specifically as follows:

[0020] ;

[0021] In the formula, , , and These correspond to the voltage vectors of the third, fifth, seventh, and ninth harmonics of the sample under test, respectively. This is the total measurement signal vector.

[0022] Preferably, the weighting factor for each harmonic is calculated as follows:

[0023] ;

[0024] In the formula, Let be the energy norm 2 of the nth harmonic system matrix. Let be the energy norm of the nth harmonic noise matrix. The full width at half maximum (FWHM) of the spread function at the nth harmonic point is... These are the weighting factors for different harmonics.

[0025] Preferably, the weighting matrix is ​​a diagonal matrix, wherein the diagonal elements are composed of weighting factors for each harmonic, specifically:

[0026] ;

[0027] In the formula, W is a diagonal matrix composed of weighting factors for different harmonics.

[0028] Preferably, the equation for solving the particle concentration distribution of the sample to be tested is:

[0029] ;

[0030] In the formula, Let v be the relaxation factor, v be the residual, and W be a diagonal matrix composed of weighting factors for different harmonics. Let c be the overall system matrix, and c be the particle concentration. This is the total measurement signal vector.

[0031] Preferably, the iterative solution of the equation includes:

[0032] The weighted algebraic reconstruction technique is used to iteratively solve the equation, where the iterative formula is:

[0033] ;

[0034] ;

[0035] In the formula, k is the number of iterations, and N is the number of columns in the system matrix. The nth element in the i-th row of the system matrix. Let i be the element of the i-th row of the voltage vector. Let i be the element of the i-th row of the residual vector. The magnetic particle concentration is obtained from the (k+1)th iteration. The particle concentration obtained from the k-th iteration is... Let i be the element of the i-th row of the voltage vector. To obtain the nth element of the particle concentration in the k-th iteration, For the i-th row of the concatenated system matrix, Let be the element of the i-th row of the residual vector, and W be a diagonal matrix composed of weighting factors for different harmonics. It is a relaxation factor.

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

[0037] (1) This invention proposes a magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting. By calculating the energy of each harmonic system matrix, the energy of the noise matrix and the FWHM of the point spread function, and combining these three parameters to form a weighting factor, the method fully considers the noise characteristics and resolution characteristics of different harmonics on the basis of traditional energy normalization weighting, so that the reconstructed image can effectively improve the spatial resolution while suppressing noise amplification.

[0038] (2) This invention fully utilizes the energy, noise, and resolution information of each harmonic by performing energy-half-peak full width joint weighted fusion of the multi-harmonic system matrix and voltage vector, thereby improving the spatial resolution of the reconstructed image and suppressing artifact generation; by normalizing the harmonic energy, the weight of higher harmonics is increased, spatial resolution is improved, and the condition number of the system matrix is ​​reduced, thus suppressing the problem of higher harmonic noise amplification in the inversion process; by normalizing the noise energy of each harmonic, the weight of harmonics affected by feedthrough and crosstalk interference is reduced, effectively suppressing artifact generation; by normalizing the full width at half peak of each harmonic, the effect of harmonics with narrower full width at half peak, i.e., higher resolution, is enhanced, effectively improving spatial resolution. Attached Figure Description

[0039] 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:

[0040] Figure 1 This is a flowchart of a magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting according to an embodiment of the present invention.

[0041] Figure 2 This is a flowchart of the weighted algebraic reconstruction algorithm according to an embodiment of the present invention. Detailed Implementation

[0042] 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.

[0043] 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.

[0044] This embodiment proposes a magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting, including:

[0045] S1. In the magnetic particle imaging system, a magnetic nanoparticle sample placed at a zero magnetic field point is excited by an excitation magnetic field, and the zero magnetic field point is moved to traverse the imaging field of view. The point spread function of several harmonics of different frequencies of the magnetic nanoparticles is measured and obtained. Based on the point spread function of each harmonic, the corresponding system matrix is ​​constructed.

[0046] S2. Under the same imaging conditions as in S1, measure a reference sample without magnetic nanoparticles and obtain noise matrices corresponding to several different frequency harmonics;

[0047] S3. Under the same imaging conditions as S1, measure the sample to be tested and obtain the magnetization response voltage vectors of several different frequency harmonics.

[0048] S4. For each harmonic, calculate the energy norm of the system matrix, the energy norm of the noise matrix, and the full width at half maximum (FWHM) of its point spread function.

[0049] S5. The harmonic system matrices obtained in S1 are spliced ​​together with the magnetization response voltage signals obtained in S3 to form the total system matrix and the total measurement signal vector; and based on the energy norm 2 of the system matrix of each harmonic, the energy norm 2 of the noise matrix, and the full width at half maximum (FWHM) of the system matrix, the weighting factor of each harmonic is calculated and a weighting matrix is ​​constructed.

[0050] S6. Using the weight matrix, the total system matrix, and the total measurement signal vector, construct the solution equation for the particle concentration distribution of the sample to be tested, and iteratively solve the solution equation to obtain the particle concentration distribution image of the sample to be tested.

[0051] Specifically, such as Figure 1 The process includes: placing the magnetic nanoparticle sample to be tested in a point-like phantom of a unit volume, moving it to the center position of the system, i.e., the zero magnetic field point, applying an excitation magnetic field to excite the magnetic nanoparticles to generate a nonlinear magnetization response, and applying a scanning field to make the zero magnetic field point traverse the entire imaging field of view; wherein, the magnetization response generated by the magnetic nanoparticles relative to different positions of the zero magnetic field point is a point spread function.

[0052] Point spread functions based on the third, fifth, seventh, and ninth harmonics of point-like magnetic nanoparticle samples. , , and Constructing a system matrix , , and ;

[0053] Pure water was placed in a point-like phantom of a unit volume, and measurements were performed using the same excitation and scanning method to obtain the noise matrices for the third, fifth, seventh, and ninth harmonics. , , and ;

[0054] The magnetic nanoparticle samples of the desired concentration were measured in the same manner to obtain the magnetization response voltage vectors for the third, fifth, seventh, and ninth harmonics. , , and ;

[0055] The system matrix of different harmonics of magnetic nanoparticles , , and and the corresponding voltage vector , , and By splicing the matrices, the overall system matrix can be obtained. and the total measurement signal vector Calculate the energy L2 norm of the system matrix for the third, fifth, seventh, and ninth harmonics respectively. , , and Energy L2 norm of the noise matrix , , and Full width at half maximum (FWHM) , , and And calculate the normalized weighting factors for different harmonics. , , and The system matrix and voltage vector obtained by splicing are reconstructed by weighted iterative calculation to obtain the particle concentration distribution c of the sample to be tested. The elements in c are correlated with the position of the zero magnetic field point to obtain the reconstructed particle concentration distribution image.

[0056] The method utilizes a self-developed magnetic nanoparticle imaging system. A displacement stage moves a unit-volume phantom containing the magnetic nanoparticle sample to the center of the system, i.e., the zero-magnetic-field point. An alternating current is then applied to the excitation coil via a power amplifier, generating an excitation magnetic field that induces a nonlinear magnetization response in the magnetic nanoparticle sample. This response is received by a receiving coil. Simultaneously, the sample is moved relative to the zero-magnetic-field point in the x-direction by the displacement stage. For each x-position, the zero-magnetic-field point is moved in the y-direction by a triangular wave scanning magnetic field, traversing the entire imaging field. The magnetization response at different positions represents the point spread function of the harmonic, yielding the point spread functions for the third, fifth, seventh, and ninth harmonics. , , and ;

[0057] Furthermore, the overall system matrix is ​​constructed by concatenating the system matrices of each harmonic row by row, specifically as follows:

[0058] ;

[0059] In the formula, The system matrix for the third harmonic of magnetic nanoparticles. The system matrix is ​​for the fifth harmonic. The system matrix for the seventh harmonic is... The system matrix is ​​for the ninth harmonic. This is the overall system matrix.

[0060] Specifically, the multi-harmonic point spread function , , and Convert to system matrix , , and Each row of the system matrix represents the magnetization response of each pixel in the entire imaging field of view when the zero magnetic field point is located at that position. Each column of the system matrix represents the magnetization response of the pixel when the zero magnetic field point traverses different positions in the imaging field of view. The conversion relationship between the system matrix and the point spread function is:

[0061] ;

[0062] In the formula, The point spread function signal of the nth harmonic of the magnetic nanoparticle at the j-th pixel when the zero magnetic field point is located at the j-th pixel; This is the system matrix.

[0063] The system matrix of different harmonics of magnetic nanoparticles , , and The overall system matrix is ​​obtained by concatenating the matrices along their rows. .

[0064] Furthermore, the different harmonic voltage vectors of the sample under test are... , , and After being converted into column vectors, they are concatenated along the row direction to obtain the total measurement signal vector. Specifically:

[0065] ;

[0066] In the formula, , , and These correspond to the voltage vectors of the third, fifth, seventh, and ninth harmonics of the sample under test, respectively. This is the total measurement signal vector.

[0067] Specifically, for a single harmonic, the voltage vector can be represented as the convolution of the point spread function and the concentration distribution of the particles to be measured. By converting the point spread function into a system matrix, the voltage vector can be represented as the product of the system matrix and the concentration distribution to be measured. For the total system matrix... and the total measurement signal vector , can be represented as:

[0068] .

[0069] Furthermore, in the field of magnetic nanoparticle imaging, low-order harmonics have high amplitudes and higher signal-to-noise ratios, but their spatial resolution is poor due to their wider point spread functions. High-order harmonics have small amplitudes and low signal-to-noise ratios, introducing numerous artifacts during reconstruction, but their point spread functions have narrower full width at half maximums (FWHMs) and richer spatial information, resulting in better spatial resolution. When performing joint reconstruction using the stitched system matrix and voltage vector, low-order harmonics, due to their larger amplitudes, account for a larger proportion in the reconstruction process, significantly limiting spatial resolution. Therefore, a redistribution of the weights of different harmonics is necessary to balance their effects.

[0070] The current mainstream weighting method is to normalize the amplitude or energy of harmonics. However, since the noise level and spatial resolution characteristics of harmonics are not considered at the same time, it is difficult to effectively improve the contribution of higher-order harmonics while suppressing noise amplification, thus failing to balance the stability of imaging and the ability to resolve details.

[0071] The method proposed in this embodiment can comprehensively consider the energy characteristics and resolution characteristics of each harmonic, effectively improving the resolution of the reconstructed image while suppressing noise amplification. Based on the stitched system matrix, the weight corresponding to each harmonic is calculated and a weight matrix is ​​constructed. Based on the weight matrix, the new system matrix, and the new measured magnetization signal voltage vector, a magnetic nanoparticle concentration distribution model is constructed. The magnetic nanoparticle concentration distribution model is iteratively reconstructed to solve for the concentration distribution of the magnetic nanoparticles to be measured.

[0072] Specifically, based on the multi-harmonic point spread function , , and The full width at half maximum (FWHM) of different harmonics was calculated. , , and The full width at half maximum (FWHM) is defined as the distance between two points in a one-dimensional point spread function where the amplitude equals the average of the maximum and minimum values.

[0073] ;

[0074] in, and x1 and x2 represent the maximum and minimum values ​​of the one-dimensional point spread function, respectively, and represent the coordinates of the two positions to the left and right of the average of the maximum and minimum values.

[0075] Noise matrices based on the third, fifth, seventh, and ninth harmonics. , , and The energy norm of the noise matrix is ​​calculated. , , and :

[0076] ;

[0077] Where m is the number of elements in the noise matrix. This represents the value of the i-th element in the nth harmonic noise matrix.

[0078] System matrix based on different harmonics , , and The energy norm of the system matrix is ​​calculated. , , and :

[0079] ;

[0080] in, Let be the value of the i-th element of the nth harmonic system matrix.

[0081] Full width at half maximum (FWHM) based on the third, fifth, seventh, and ninth harmonics. , , and Energy L2 norm of the noise matrix , , and and the energy norm of the system matrix , , and The weighting factors for different harmonics are obtained:

[0082] ;

[0083] In the formula, Let be the energy norm 2 of the nth harmonic system matrix. Let be the energy norm of the nth harmonic noise matrix. The full width at half maximum (FWHM) of the spread function at the nth harmonic point is... These are the weighting factors for different harmonics.

[0084] Furthermore, the weight matrix is ​​a diagonal matrix, where the diagonal elements are composed of weighting factors for each harmonic, specifically:

[0085] ;

[0086] In the formula, W is a diagonal matrix composed of weighting factors for different harmonics.

[0087] Furthermore, based on the obtained overall system matrix Total measurement signal vector And with the weight matrix W, the equation for solving the particle concentration distribution of the sample to be tested is as follows:

[0088] ;

[0089] In the formula, Let v be the relaxation factor, v be the residual, and W be a diagonal matrix composed of weighting factors for different harmonics. Let c be the overall system matrix, and c be the particle concentration. This is the total measurement signal vector.

[0090] The Algebraic Reconstruction Algorithm (ART) is used to iteratively reconstruct the concentration distribution of the magnetic nanoparticles to be solved. The iterative formula is as follows:

[0091] ;

[0092] ;

[0093] In the formula, k is the number of iterations, and N is the number of columns in the system matrix. The nth element in the i-th row of the system matrix. Let i be the element of the i-th row of the voltage vector. Let i be the element of the i-th row of the residual vector. The magnetic particle concentration is obtained from the (k+1)th iteration. The particle concentration obtained from the k-th iteration is... Let i be the element of the i-th row of the voltage vector. To obtain the nth element of the particle concentration in the k-th iteration, For the i-th row of the concatenated system matrix, Let be the element of the i-th row of the residual vector, and W be a diagonal matrix composed of weighting factors for different harmonics. It is a relaxation factor.

[0094] This embodiment introduces an energy-full width at half maximum (FWHM) joint weighting mechanism that simultaneously considers harmonic energy characteristics and spatial resolution characteristics. This mechanism leverages the high signal-to-noise ratio of low-order harmonics while fully utilizing the high resolution of high-order harmonics, effectively improving the stability of the system matrix reconstruction and the fidelity of image details without amplifying noise artifacts. This method overcomes the limitations of traditional weighting strategies based solely on energy normalization, achieving efficient utilization of multi-harmonic information, significantly improving the spatial resolution of reconstructed images, and suppressing artifacts. This provides reliable technical support for the widespread application of magnetic nanoparticle imaging in precision medicine, early lesion detection, and cell tracking.

[0095] To more clearly illustrate the technical solution of the present invention, specific embodiments are provided below for description:

[0096] The practical application of this invention is an extension of single-harmonic narrowband magnetic nanoparticle imaging, and this will be explained using single-harmonic narrowband magnetic nanoparticle imaging as an example.

[0097] First, the magnetic nanoparticle sample to be tested is loaded into a unit volume point-like biomimetic force and moved to the center of the imaging system by a displacement stage. The excitation magnetic field strength and imaging field width are set. These parameters determine the magnitude of the excitation magnetic field and the scanning magnetic field. After the zero magnetic field point traverses the entire imaging field, the point spread functions of the third, fifth, seventh and ninth harmonics of the magnetization response signal are extracted.

[0098] Pure water was placed in a unit volume point phantom and moved to the center of the system. Measurements were performed using the same excitation and scanning parameters to obtain noise matrices for different harmonics.

[0099] The magnetic nanoparticle sample to be measured is placed in the imaging field of view. The sample is excited under the same excitation conditions, and the voltage vectors of the third, fifth, seventh, and ninth harmonics of its magnetization response signal are extracted. Taking the third harmonic as an example, the PSF and the voltage vector can be represented by convolution:

[0100] ;

[0101] In the formula, This is the voltage vector of the third harmonic magnetization response. Let be the point spread function of the third harmonic, and c be the particle concentration distribution.

[0102] By using FFP point meshing, the convolution relationship can be transformed into a matrix multiplication relationship. Taking the third harmonic as an example, the formula for transforming PSF into the system matrix S is:

[0103] ;

[0104] in, It is the magnetization response of the third harmonic of the i-th pixel in the imaging field of view when the position of FFP is the j-th pixel, and M is the number of imaging pixels.

[0105] Taking the third harmonic as an example, the system matrix and voltage vector The relationship between them can be represented as:

[0106] ;

[0107] Similarly, the same expression can be used for the fifth, seventh, and ninth harmonics, where 'c' is the same, representing the particle concentration distribution to be solved.

[0108] By fusing and jointly inverting the system matrices and voltage vectors of low-order and high-order harmonics, the advantages of high signal-to-noise ratio of low-order harmonics and high resolution of high-order harmonics can be combined. The multi-harmonic system matrix fusion equation can be expressed as:

[0109] ;

[0110] in, The system matrix for the third harmonic of magnetic nanoparticles. The system matrix is ​​for the fifth harmonic. The system matrix for the seventh harmonic is... This is the system matrix for the ninth harmonic; note that the elements in the system matrix contain both real and imaginary parts. , , and These correspond to the voltage vectors of the third, fifth, seventh, and ninth harmonics of the sample under test, respectively. The fused system matrix and voltage vectors can be expressed as follows: and :

[0111] , ;

[0112] In practical systems, the intensity of lower harmonics is greater than that of higher harmonics; for example, the intensity of the third harmonic is approximately ten times that of the ninth harmonic. Weighting factors need to be assigned to the system matrix for different harmonics to balance their interactions. Compared to traditional energy-normalized weighting, this embodiment proposes an energy-full width at half maximum (FWHM) joint weighting, which fully considers the energy, noise, and resolution characteristics of different harmonics, and the weighting factors for different harmonics... :

[0113] ;

[0114] in, Let be the energy norm 2 of the nth harmonic system matrix. Let be the energy norm of the nth harmonic noise matrix. Let be the full width at half maximum (FWHM) of the spread function at the nth harmonic point.

[0115] By calculating the weighting factors for different harmonics, the weighting matrix can be obtained. :

[0116] ;

[0117] Where W is a diagonal matrix composed of weighting factors for different harmonics.

[0118] The resulting spliced ​​system matrix Voltage vector With the weight matrix W, a model can be constructed to solve for the concentration distribution of the magnetic nanoparticles to be measured:

[0119] ;

[0120] in, is the relaxation factor, and v is the residual.

[0121] The weighted algebraic reconstruction algorithm (ART) can be used to iteratively reconstruct the concentration distribution of magnetic nanoparticles to be solved. The specific algorithm flow is as follows: Figure 2 As shown:

[0122] First, input the point spread functions for different harmonics. , , and noise matrix , , and and voltage vector , , and Input the maximum number of iterations k, and the convergence error. and relaxation factor ;

[0123] Transform the multiharmonic point spread function into a system matrix. , , and The system matrix is ​​formed by concatenating the corresponding voltage vectors. and the total measurement signal vector And calculate the energy norm of the system matrix for different harmonics. , , and Energy L2 norm of the noise matrix , , and and half-peak full width , , and Obtain the weighting factors for different harmonics. And the weight matrix W.

[0124] Initialize the residual v=0, and perform traversal projection based on the number of system equations. If the concentration distribution c obtained after a certain traversal satisfies the convergence error, the iteration process is terminated; otherwise, iteration continues until the maximum number of iterations is reached. The specific iteration formula is as follows:

[0125] ;

[0126] ;

[0127] Where k is the number of iterations, and N is the number of columns in the system matrix. The nth element in the i-th row of the system matrix. Let i be the element of the i-th row of the voltage vector. Let i be the element of the i-th row of the residual vector.

[0128] This embodiment proposes a magnetic particle imaging method based on multi-harmonic energy-full-width-of-half (HWHM) joint weighting. By jointly inverting multiple harmonics, it combines the advantages of high signal-to-noise ratio of low-order harmonics and high resolution of high-order harmonics, resulting in reconstructed images with good spatial resolution and fewer artifacts. Compared with traditional energy normalization weighting, the energy-full-width-of-half (HWHM) weighting fully utilizes the noise and resolution characteristics of harmonics, significantly improving spatial resolution without introducing excessive noise. Moreover, the weight matrix reduces the condition number and improves the orthogonality of the system matrix, effectively reducing the amplification effect of high-order harmonic noise during inverse problem solving and accelerating the convergence speed of iterative reconstruction. This method is of great significance for advancing the clinical application of magnetic nanoparticle imaging.

[0129] The above are merely preferred embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting, characterized in that, include: S1. In the magnetic particle imaging system, a magnetic nanoparticle sample placed at a zero magnetic field point is excited by an excitation magnetic field, and the zero magnetic field point is moved to traverse the imaging field of view. The point spread function of several harmonics of different frequencies of the magnetic nanoparticles is measured and obtained. Based on the point spread function of each harmonic, the corresponding system matrix is ​​constructed. S2. Under the same imaging conditions as in S1, measure a reference sample without magnetic nanoparticles and obtain noise matrices corresponding to several different frequency harmonics; S3. Under the same imaging conditions as S1, measure the sample to be tested and obtain the magnetization response voltage vectors of several different frequency harmonics. S4. For each harmonic, calculate the energy norm of the system matrix, the energy norm of the noise matrix, and the full width at half maximum (FWHM) of its point spread function. S5. The harmonic system matrices obtained in S1 are spliced ​​together with the magnetization response voltage signals obtained in S3 to form the total system matrix and the total measurement signal vector; and based on the energy norm 2 of the system matrix of each harmonic, the energy norm 2 of the noise matrix, and the full width at half maximum (FWHM) of the system matrix, the weighting factor of each harmonic is calculated and a weighting matrix is ​​constructed. S6. Using the weight matrix, the total system matrix, and the total measurement signal vector, construct the solution equation for the particle concentration distribution of the sample to be tested, and iteratively solve the solution equation to obtain the particle concentration distribution image of the sample to be tested.

2. The magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting according to claim 1, characterized in that, The point spread functions of the magnetic nanoparticles for several different frequency harmonics include the point spread functions of the third, fifth, seventh, and ninth harmonics; the noise matrices corresponding to the several different frequency harmonics include the noise matrices of the third, fifth, seventh, and ninth harmonics.

3. The magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting according to claim 2, characterized in that, The overall system matrix is ​​constructed by concatenating the system matrices of each harmonic row by row, specifically as follows: ； In the formula, The system matrix for the third harmonic of magnetic nanoparticles. The system matrix is ​​for the fifth harmonic. The system matrix for the seventh harmonic is... The system matrix is ​​for the ninth harmonic. This is the overall system matrix; The total measurement signal vector is constructed by concatenating the magnetization response voltage signal vectors of each harmonic row by row, specifically as follows: ； In the formula, , , and These correspond to the voltage vectors of the third, fifth, seventh, and ninth harmonics of the sample under test, respectively. This is the total measurement signal vector.

4. The magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting according to claim 1, characterized in that, The weighting factors for each harmonic are calculated as follows: ； In the formula, Let be the energy norm 2 of the nth harmonic system matrix. Let be the energy norm of the nth harmonic noise matrix. The full width at half maximum (FWHM) of the spread function at the nth harmonic point is... These are the weighting factors for different harmonics.

5. The magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting according to claim 4, characterized in that, The weight matrix is ​​a diagonal matrix, where the diagonal elements are composed of weighting factors for each harmonic, specifically: ； In the formula, W is a diagonal matrix composed of weighting factors for different harmonics.

6. The magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting according to claim 1, characterized in that, The equation for solving the particle concentration distribution of the sample to be tested is: ； In the formula, Let v be the relaxation factor, v be the residual, and W be a diagonal matrix composed of weighting factors for different harmonics. Let c be the overall system matrix, and c be the particle concentration. This is the total measurement signal vector.

7. The magnetic particle imaging method based on multi-harmonic energy-full width at half maximum (FWHM) joint weighting according to claim 1, characterized in that, The equation is solved iteratively, including: The weighted algebraic reconstruction technique is used to iteratively solve the equation, where the iterative formula is: ； ； In the formula, k is the number of iterations, and N is the number of columns in the system matrix. The nth element in the i-th row of the system matrix. Let i be the element of the i-th row of the voltage vector. Let i be the element of the i-th row of the residual vector. The magnetic particle concentration is obtained from the (k+1)th iteration. The particle concentration obtained from the k-th iteration is... Let i be the element of the i-th row of the voltage vector. To obtain the nth element of the particle concentration in the k-th iteration, For the i-th row of the concatenated system matrix, Let be the element of the i-th row of the residual vector, and W be a diagonal matrix composed of weighting factors for different harmonics. It is a relaxation factor.

Images