A method and system for constructing a lesion model based on a magnetocardiogram

Abstract

The application discloses a kind of based on magnetocardiogram's lesion model construction method and system, comprising: acquisition multi-channel magnetocardiogram and preprocessing are obtained preprocessing observation matrix;Based on second-order statistics blind source separation of second-order statistics blind source separation based on multi-delay covariance weighting joint diagonalization, generate interference suppression weight and reconstruction are obtained purification observation matrix;Coupling separation matrix and interference suppression weight generate inversion prior parameter and establish weighted forward model;Under the constraint of weighted forward model, using the improved FOCUSS of connected domain focus to the sparse reconstruction of purification observation matrix is obtained convergent source distribution;Based on convergent source distribution extraction lesion voxel boundary and reconstruction surface grid, combined with energy projection consistency check output three-dimensional lesion model.The application improves lesion positioning accuracy and three-dimensional model construction stability under strong noise and individual difference conditions.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional modeling technology, and in particular to a method and system for constructing lesion models based on magnetocardiography. Background Technology

[0002] Magnetocardiography (MCG) is used to non-invasively record the magnetic field distribution of cardiac electrical activity in the extracellular space. It features high temporal resolution and relatively minimal impact on tissue conductivity, and has been used for arrhythmia localization, ischemia assessment, and preoperative planning. Current techniques typically use multi-channel observation signals as input, first performing bandpass filtering, power frequency suppression, and baseline correction. Then, a forward mapping model between the cardiac source and the observation is established by combining sensor geometric information. Under the constraints of this model, an underdetermined inverse problem is solved to reconstruct the source distribution.

[0003] The inversion process includes minimum norm inversion, beamforming, and dipole fitting methods. Minimum norm inversion obtains a smooth source distribution by minimizing energy, but is prone to diffusion and multi-peak interference. Beamforming relies on spatial filters to suppress non-target directional components, and is susceptible to sidelobe artifacts when affected by noise and model mismatch. Dipole fitting depends on initial values ​​and model assumptions, and lacks stability when dealing with multifocal or complex lesions. To improve focality, some schemes introduce sparse constraints and iterative reweighting solutions; however, under conditions of strong magnetocardiographic noise, significant individual variability, and unstable statistical structure of the observation matrix, sparse inversion is easily influenced by artifacts, leading to false localization and non-reproducible results.

[0004] In signal purification and 3D model construction, existing solutions often employ blind source separation based on independence or second-order statistics to suppress respiratory electromyography, environmental magnetic noise, and sensor drift. The purified observations are then input into sparse inversion, and the lesion model is obtained through threshold segmentation, connected component aggregation, and mesh reconstruction. However, due to the lack of coupling constraints between blind source separation and sparse inversion for inversion stability, and insufficient verification of mesh and source distribution consistency, model boundaries are sensitive to thresholds and noise. Holes and leaking voxels easily lead to surface mesh distortion and positioning errors. Furthermore, the relative pose of the magnetocardiogram sensor array and the anatomical location of the heart varies, and the determination of source space boundaries and voxel modeling depends on registration information. Registration errors amplify inversion uncertainties and affect subsequent boundary extraction.

[0005] Therefore, how to provide a method and system for constructing lesion models based on magnetocardiography is a problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0006] One objective of this invention is to propose a method and system for constructing lesion models based on magnetocardiography. This invention achieves stable lesion localization and outputs a three-dimensional model with consistency verification by using blind source separation and purification and connected domain focusing sparse inversion coupled modeling.

[0007] A method for constructing a lesion model based on magnetocardiography according to an embodiment of the present invention includes:

[0008] Raw observation data from multi-channel magnetocardiography were collected and preprocessed to obtain a preprocessed observation matrix;

[0009] The set of multi-delay covariance matrices is calculated based on the preprocessed observation matrix, and joint diagonalization is performed by the second-order statistical blind source separation algorithm to obtain the separation matrix and source component sequence. Interference suppression weights are generated based on the autocorrelation peak and spatial projection energy ratio of the source component sequence to obtain the purified observation matrix.

[0010] The separation matrix and the interference suppression weights are coupled to generate inversion prior parameters, and a weighted forward model is established based on the inversion prior parameters.

[0011] Under the constraints of the weighted forward model, the improved FOCUSS algorithm is used to perform sparse reconstruction processing on the cleaned observation matrix to obtain the convergence source distribution;

[0012] The convergence source distribution is localized, and a set of candidate lesion voxels is generated based on the peak connected domain of the convergence source distribution. The boundary of the candidate lesion voxel set is determined based on the energy gradient of the convergence source distribution, and the lesion voxel boundary set is obtained.

[0013] Three-dimensional voxel aggregation processing is performed on the voxel boundary set of the lesion, morphological closing operation is used to eliminate holes and surface triangulation reconstruction is performed to obtain the three-dimensional lesion surface mesh;

[0014] Perform consistency verification on the three-dimensional lesion surface mesh and the convergence source distribution, calculate the energy projection value of the source distribution corresponding to the mesh vertex and remove mesh vertices below the threshold, and output the three-dimensional lesion model.

[0015] Optionally, the preprocessing to obtain a preprocessed observation matrix is ​​as follows:

[0016] Collect raw observation data of multi-channel magnetocardiography and simultaneously record the spatial coordinates of each channel sensor, sampling rate parameters and channel number information to obtain the raw observation matrix and channel coordinate table;

[0017] Based on the channel coordinate table, a spatial basis function is constructed. For each sampling time of the original observation matrix, the spatial basis function is fitted on the entire channel. The fitting result is then subtracted from the original observation matrix to obtain the far-field observation matrix.

[0018] In the far-field observation matrix, the band-limited coherence index between each channel and its neighboring channel set is calculated. Based on the band-limited coherence index, a channel reliability weight vector is generated and applied to the far-field observation matrix to obtain the reliability-weighted observation matrix.

[0019] The heartbeat marker sequence is calculated based on the reliability-weighted observation matrix, and the reliability-weighted observation matrix is ​​aligned by periodic slicing and periodic resampling to obtain the aligned periodic observation matrix.

[0020] The sample confidence mask matrix is ​​calculated based on the aligned periodic observation matrix and applied to the aligned periodic observation matrix to obtain the robust periodic observation matrix;

[0021] The periodic dimension weighted aggregation process is performed on the robust periodic observation matrix to obtain the continuous robust observation matrix. The noise covariance matrix is ​​calculated for the sample intervals with energy envelopes below the threshold. The continuous robust observation matrix is ​​then whitened based on the noise covariance matrix to obtain the preprocessed observation matrix.

[0022] Optionally, obtaining the purified observation matrix specifically involves:

[0023] The cross-channel energy envelope sequence is calculated based on the preprocessed observation matrix, and autocorrelation analysis is performed to obtain the estimated value of the dominant cardiac cycle and the corresponding delay value range.

[0024] Multiple delay parameters are selected within the delay value range. The covariance matrix of the corresponding delay is calculated for the preprocessed observation matrix. Matrix confidence weights are generated based on the diagonal dominance index of each covariance matrix to obtain a set of multi-delay covariance matrices and a set of matrix confidence weights.

[0025] The set of multi-delay covariance matrices and the set of matrix confidence weights are used as inputs. The second-order statistical blind source separation algorithm is used to perform weighted joint diagonalization processing to output the separation matrix. The source component sequence is obtained by linear transformation of the preprocessed observation matrix by the separation matrix. At the same time, the equivalent mixing matrix consistent with the separation matrix is ​​calculated.

[0026] For each source component sequence, calculate the set of autocorrelation peak indices within the time delay range;

[0027] Based on the column vectors of the equivalent mixing matrix, a set of basis vectors for the source space is constructed. The projection energy of each column vector on the basis vector set for the source space and the residual energy on the orthogonal complement space are calculated respectively to obtain a set of spatial projection energy ratio indices.

[0028] Monotonic normalization is performed on the set of autocorrelation peak indices and the set of spatial projection energy ratio indices to generate interference suppression weights;

[0029] The source component sequence is subjected to component-level weighted suppression based on the interference suppression weight, and the weighted source component sequence is reconstructed through an equivalent mixing matrix to output a cleaned observation matrix.

[0030] Optionally, the coupling of the separation matrix and the interference suppression weights to generate inversion prior parameters, and the establishment of a weighted forward model based on the inversion prior parameters, specifically involves:

[0031] The source component sequences are sorted according to the interference suppression weights to generate a set of valid source component indices;

[0032] Based on the set of valid source component indices, the corresponding column vectors are extracted from the separation matrix to form a valid separation submatrix. The column vectors of the valid separation submatrix are then scaled to obtain the coupling separation matrix.

[0033] The reliability vector of each channel is calculated based on the coupling separation matrix, and the interference suppression weights corresponding to the effective source component index set are monotonically normalized to generate the component confidence sequence.

[0034] The channel reliability vector and component confidence sequence are combined to generate inversion prior parameters, and the channel weighting factor is generated based on the channel reliability vector, and the source space weighting factor is generated based on the component confidence sequence.

[0035] The spatial coordinates of each channel sensor recorded during multi-channel magnetocardiography acquisition are obtained, and a source space voxel coordinate set is generated within the preset cardiac source space boundary. The mapping coefficient of each voxel to the observation contribution of each channel sensor is calculated based on the spatial coordinates of each channel sensor and the source space voxel coordinate set, forming an unweighted forward model.

[0036] The channel dimension of the unweighted forward model is weighted according to the channel weighting factor, and the voxel dimension of the unweighted forward model is weighted according to the source space weighting factor. Then, the column vector normalization process is performed on the coefficient vector corresponding to each voxel after weighting to obtain the weighted forward model.

[0037] Optionally, the improved FOCUSS algorithm is specifically as follows:

[0038] Under the weighted forward model constraints, the minimum norm solution of the cleaned observation matrix that satisfies the observation fitting constraints is obtained to generate the initial source distribution;

[0039] Calculate the voxel energy corresponding to each voxel in the source space voxel coordinate set based on the initial source distribution, and establish a three-dimensional adjacency relationship on the source space voxel coordinate set to generate a voxel adjacency table.

[0040] Select the set of voxels with energy peaks based on voxel energy, and use the set of voxels with energy peaks as the starting point for region growth. Generate a set of connected candidate regions under the constraint of voxel adjacency list according to the expansion order of voxel energy from high to low.

[0041] For each connected candidate domain in the set of connected candidate domains, calculate the domain energy and the number of domain voxels, and calculate the domain-level confidence weights based on the domain energy and the number of domain voxels to obtain the set of domain-level confidence weights;

[0042] A connected-domain focused sparse weight vector is generated based on the domain-level confidence weight set, and the voxel dimension of the weighted forward model is scaled to obtain the connected-domain weighted forward model.

[0043] Under the constraints of the connected component weighted forward model, the weighted least squares subproblem of cleaning the observation matrix is ​​solved to obtain the intermediate solution vector. The intermediate solution vector is then merged and updated with the connected component focused sparse weight vector to generate the current source distribution.

[0044] The residual matrix is ​​calculated based on the weighted forward model of connected components and the current source distribution. The ratio of the number of elements in the intersection to the number of elements in the union of the candidate connected component sets in two adjacent rounds is also calculated. The energy change of the residual matrix in two adjacent rounds is calculated to generate the convergence determination quantity.

[0045] When the convergence determination quantity meets the preset termination threshold, the convergence source distribution is output. When the convergence determination quantity does not meet the preset termination threshold, the current source distribution is used as the initial source distribution for the next round and the iteration count is updated. When the maximum number of iterations is reached, the convergence source distribution is output.

[0046] Optionally, the localization process performed on the convergence source distribution specifically includes:

[0047] Perform voxel energy normalization on the convergence source distribution, calculate the voxel energy statistics of the convergence source distribution, and generate energy quantile threshold parameters and normalized energy voxel maps.

[0048] Voxels whose voxel energies are not lower than the energy quantile threshold parameter in the normalized energy voxel graph are marked as high-energy voxels. A three-dimensional six-adjacency relationship is established on the source space voxel coordinate set, and a set of peak connected regions is generated based on the high-energy voxels.

[0049] Calculate the domain confidence weight set based on the peak connected component set;

[0050] Based on the set of domain confidence weights, peak connected regions with domain confidence weights greater than a preset threshold are selected. The selected peak connected regions are merged to obtain a set of candidate lesion voxels. A layer of adjacent voxels is then extended outside the set of candidate lesion voxels to generate a boundary search voxel band.

[0051] The spatial energy gradient vector field of the convergent source distribution is calculated on the boundary search voxel band, and non-maximum suppression is performed to obtain the candidate boundary voxel set.

[0052] The candidate boundary voxel set is subjected to closed connectivity verification, boundary segments that do not form closed toruses are deleted, and the gap voxels are subjected to connectivity completion processing to obtain the lesion voxel boundary set.

[0053] Optionally, the three-dimensional voxel aggregation process performed on the lesion voxel boundary set specifically includes:

[0054] Map the set of voxel boundaries of the lesion to the set of voxel coordinates of the source space to generate a boundary occupation marker voxel map and a set of boundary voxels.

[0055] Based on the three-dimensional six-adjacency relationship of the source space voxel coordinate set, the boundary voxel set is subjected to connectivity aggregation processing, and the parity counting and internal / external discrimination are performed along three orthogonal directions to obtain the lesion entity voxel set;

[0056] Calculate the discrete shortest distance from the voxel to the boundary occupied marker voxel map based on the voxel set of lesion entities, generate a discrete distance voxel map, and extract the candidate voxel set of holes from the discrete distance voxel map.

[0057] Based on the discrete distance voxel map, structural element scale parameters are generated at the corresponding positions of the candidate voxel set of pores, and three-dimensional morphological closing operations are performed on the lesion entity voxel set to obtain the pore repair voxel set.

[0058] The sign distance value of the voxel corner points is calculated based on the voxel set of hole repair, and the intersection point is determined at the voxel boundary to obtain the initial three-dimensional lesion surface mesh;

[0059] The non-manifold edges and boundary loops of the initial 3D lesion surface mesh are detected, and boundary loop filling is performed. Local Laplacian smoothing is then applied to the filled area to output the 3D lesion surface mesh.

[0060] Optionally, the consistency verification process for the three-dimensional lesion surface mesh and the convergence source distribution is specifically as follows:

[0061] The coordinates of the grid vertices of the three-dimensional lesion surface mesh are transformed to the source space voxel coordinate set coordinate system. For each grid vertex, the voxel unit containing the grid vertex is determined, and a neighborhood voxel set and a vertex sampling index table are generated.

[0062] The voxel energy of the neighborhood voxel set in the convergence source distribution is read from the vertex sampling index table, and a weighted summation is performed according to the distance attenuation weight from the voxel center to the grid vertex to generate an energy projection value sequence.

[0063] Calculate the projection value statistics and projection value quantile threshold parameters based on the energy projection value sequence. Mark the grid vertices whose projection value statistics are lower than the projection value quantile threshold parameters as low-consistency vertices to obtain a set of low-consistency vertices.

[0064] In the 3D lesion surface mesh, delete the triangles adjacent to the set of low-consistency vertices while maintaining the connectivity of the remaining triangles to generate a trimmed lesion surface mesh.

[0065] Based on the mesh of the lesion surface, a mesh-occupied voxel map is constructed, and voxels whose voxel energy is not lower than the energy quantile threshold parameter in the convergence source distribution are marked as high-energy voxel sets. The proportion of voxels falling outside the mesh-occupied voxel map of the high-energy voxel sets is calculated to generate an external leakage consistency index.

[0066] When the leakage consistency index is higher than the preset leakage threshold, delete the triangular pieces adjacent to the leaked voxels outside the grid on the voxel map, and perform boundary ring filling to obtain a three-dimensional lesion model.

[0067] Optionally, a lesion model construction system based on magnetocardiography includes:

[0068] The multi-channel acquisition and coordinate recording module acquires raw observation data of multi-channel magnetocardiogram and records sensor spatial coordinates, sampling rate parameters and channel number information, and outputs raw observation matrix and channel coordinate table;

[0069] The preprocessing module performs spatial basis function fitting to subtract far-field interference, band-limited coherence weighting, heart rate cycle alignment, and noise covariance whitening on the original observation matrix, and outputs the preprocessed observation matrix.

[0070] The weighted joint diagonalization purification module constructs a set of multi-delay covariance matrices and a set of matrix confidence weights, performs weighted joint diagonalization processing, and outputs a purified observation matrix.

[0071] The inversion prior coupling and weighted forward modeling module performs coupling on the separation matrix and interference suppression weights to generate inversion prior parameters, generates channel weighting factors and source space weighting factors, and establishes a weighted forward model;

[0072] The connected domain focusing FOCUSS sparse reconstruction module generates a set of connected candidate domains based on voxel energy and three-dimensional adjacency relations, calculates domain-level confidence weights, constructs a connected domain weighted forward model, and iteratively solves for the output convergence source distribution.

[0073] The focalization and boundary extraction module generates a set of peak connected regions and a set of candidate lesion voxels based on the energy quantile threshold parameter, calculates the spatial energy gradient vector field and performs non-maximum suppression and closed connectivity verification, and outputs a set of lesion voxel boundaries.

[0074] The mesh reconstruction and consistency verification module performs three-dimensional voxel aggregation, morphological closing operation and surface triangulation reconstruction on the lesion voxel boundary set to obtain a three-dimensional lesion surface mesh, calculates the energy projection value sequence and leakage consistency index and trims and repairs it, and outputs a three-dimensional lesion model.

[0075] The beneficial effects of this invention are:

[0076] (1) This invention achieves statistical structure stabilization and channel consistency enhancement of the preprocessed observation matrix by performing spatial basis function fitting to deduct far-field interference, band-limited coherence weighting, heart rate cycle alignment and noise covariance whitening on the original observation matrix. This effectively improves the usability of constructing multi-delay covariance matrix sets and weighted joint diagonalization processing, and improves the signal-to-noise ratio and repeatability of the purified observation matrix under strong noise and individual differences.

[0077] (2) This invention constructs a set of multi-delay covariance matrices and introduces a set of matrix confidence weights to perform weighted joint diagonalization. At the same time, it generates interference suppression weights based on the set of autocorrelation peak index and the set of spatial projection energy ratio index and reconstructs the purified observation matrix. This realizes the distinction and suppression of cardiac source-related components and interference components, effectively reduces the probability of artifact components entering the inversion link, and enhances the stability and consistency of sparse reconstruction results under the constraints of subsequent weighted forward model.

[0078] (3) The improved FOCUSS algorithm in this invention generates a connected domain focused sparse weight vector and constructs a connected domain weighted forward model by introducing a voxel energy-driven connected candidate domain set and a domain-level confidence weight set. This achieves focused convergence of the sparse support set to the connected focal region, significantly suppresses discrete spurious lesions and improves the spatial coherence of the convergence source distribution and the lesion voxel boundary set. As a result, a more stable three-dimensional lesion model output is obtained in the three-dimensional voxel aggregation, surface triangulation reconstruction and consistency verification processing. Attached Figure Description

[0079] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings:

[0080] Figure 1 This is a flowchart of a method for constructing a lesion model based on magnetocardiography proposed in this invention;

[0081] Figure 2 This is a system structure diagram of a lesion model construction system based on magnetocardiography proposed in this invention;

[0082] Figure 3 This is a data flow diagram of the improved FOCUSS algorithm for a lesion model construction method based on magnetocardiography proposed in this invention. Detailed Implementation

[0083] Example 1: Reference Figures 1-3 A method for constructing a lesion model based on magnetocardiography, comprising:

[0084] Raw observation data from multi-channel magnetocardiography were collected and preprocessed to obtain a preprocessed observation matrix;

[0085] The set of multi-delay covariance matrices is calculated based on the preprocessed observation matrix, and joint diagonalization is performed by the second-order statistical blind source separation algorithm to obtain the separation matrix and source component sequence. Interference suppression weights are generated based on the autocorrelation peak and spatial projection energy ratio of the source component sequence to obtain the purified observation matrix.

[0086] The separation matrix and the interference suppression weights are coupled to generate inversion prior parameters, and a weighted forward model is established based on the inversion prior parameters.

[0087] Under the constraints of the weighted forward model, the improved FOCUSS algorithm is used to perform sparse reconstruction processing on the cleaned observation matrix to obtain the convergence source distribution;

[0088] The convergence source distribution is localized, and a set of candidate lesion voxels is generated based on the peak connected domain of the convergence source distribution. The boundary of the candidate lesion voxel set is determined based on the energy gradient of the convergence source distribution, and the lesion voxel boundary set is obtained.

[0089] Three-dimensional voxel aggregation processing is performed on the voxel boundary set of the lesion, morphological closing operation is used to eliminate holes and surface triangulation reconstruction is performed to obtain the three-dimensional lesion surface mesh;

[0090] Perform consistency verification on the three-dimensional lesion surface mesh and the convergence source distribution, calculate the energy projection value of the source distribution corresponding to the mesh vertex and remove mesh vertices below the threshold, and output the three-dimensional lesion model.

[0091] In this embodiment, the preprocessing to obtain the preprocessed observation matrix specifically involves:

[0092] Collect raw observation data of multi-channel magnetocardiography and simultaneously record the spatial coordinates of each channel sensor, sampling rate parameters and channel number information to obtain the raw observation matrix and channel coordinate table;

[0093] Based on the channel coordinate table, a spatial basis function is constructed. For each sampling time of the original observation matrix, the spatial basis function is fitted on the entire channel. The fitting result is then subtracted from the original observation matrix to obtain the far-field observation matrix.

[0094] In this embodiment 1, the channel coordinate table records the three coordinate components of the sensor's spatial coordinates row by row according to the channel number; the mean of all channels is calculated for each of the three coordinate components in the channel coordinate table, and the mean of the corresponding coordinate component of each channel is subtracted to complete the coordinate centering; the Euclidean distance to the origin is calculated based on the centered coordinates of each channel, and the maximum distance of all channels is taken. The maximum distance of all channels is used to normalize the three centered coordinate components; the spatial basis function is defined on the normalized coordinates as 1 constant term, 3 linear coordinate terms, 3 coordinate self-multiplication terms, and 3 coordinate pairwise product terms; the normalized coordinates of each channel are substituted into the spatial basis function and 10 basis function values ​​are calculated and arranged in order of channel number to form the spatial basis function;

[0095] In the far-field observation matrix, the band-limited coherence index between each channel and its neighboring channel set is calculated. Based on the band-limited coherence index, a channel reliability weight vector is generated and applied to the far-field observation matrix to obtain the reliability-weighted observation matrix.

[0096] In this embodiment 1, the Euclidean distance between the target channel and other channel sensor spatial coordinates is calculated based on the channel coordinate table, and the six channels with the smallest distance are selected to form a neighborhood channel set. For the target channel and the corresponding channel in the neighborhood channel set in the far-field observation matrix, a piecewise Fourier transform with a window length equal to the number of samples per second corresponding to the sampling rate is used, and the energy is averaged on the piecewise results to obtain the autospectrum of each channel and the cross-spectrum between the target channel and the neighborhood channels. The average spectrum is obtained by averaging the autospectrum of all channels and locating the maximum peak frequency of the average spectrum. The lower and upper band-limit boundaries are determined by searching along both sides of the peak frequency point until the spectral amplitude drops to 0.5 of the peak spectral amplitude. The ratio of the cross-spectral energy to the product of the autospectral energy of the two channels is calculated within the frequency range defined by the lower and upper band-limit boundaries and averaged in the frequency dimension to obtain the band-limited coherence index of the target channel relative to each neighborhood channel. The arithmetic mean of the band-limited coherence indices of the target channel and the neighborhood channel set is taken to obtain the channel reliability index. The minimum and maximum values ​​of all channel reliability indices are calculated and linearly normalized to obtain a channel reliability weight vector from 0 to 1.

[0097] The heartbeat marker sequence is calculated based on the reliability-weighted observation matrix, and the reliability-weighted observation matrix is ​​aligned by periodic slicing and periodic resampling to obtain the aligned periodic observation matrix.

[0098] In this embodiment 1, the energy of the entire channel is summarized point by point on the reliability-weighted observation matrix to form a global energy sequence; the global energy sequence is subjected to mean removal processing and the autocorrelation sequence is calculated. The sample interval corresponding to the first main peak of the autocorrelation sequence after zero delay is located and used as the cardiac cycle sample interval; a moving average smoother is constructed using half of the cardiac cycle sample interval and the global energy sequence is smoothed to obtain a smoothed energy sequence; the difference between adjacent samples is calculated on the smoothed energy sequence and the inflection point where the sign of the difference changes from positive to negative is located as a candidate peak point; an amplitude discriminant is constructed based on the median and the absolute deviation of the median of the smoothed energy sequence, and candidate peak points lower than the amplitude discriminant are eliminated to obtain a set of effective peak points; the set of effective peak points is sorted by time, and peak points with larger amplitudes are retained when the interval between adjacent peak points is less than half of the cardiac cycle sample interval, and the peak point sampling sequence number is output to form a cardiac marker sequence;

[0099] The sample confidence mask matrix is ​​calculated based on the aligned periodic observation matrix and applied to the aligned periodic observation matrix to obtain the robust periodic observation matrix;

[0100] In this embodiment 1, the aligned period observation matrix is ​​organized according to the cardiac cycle index, channel index, and aligned sampling number, while maintaining consistency between the period indexes. For each channel index and each aligned sampling number, all cardiac cycle samples are summarized, and the median and median absolute deviation are calculated. For each cardiac cycle sample, the absolute value of the difference from the median is calculated, and the median absolute deviation is used to complete the normalization to obtain the deviation metric. Within the same channel index and the same aligned sampling number, all deviation metrics are sorted, and the 90th percentile deviation metric is taken as the cut-off point. Cardiac cycle samples with a deviation metric not greater than the cut-off point are assigned a value of 1, and cardiac cycle samples with a deviation metric greater than the cut-off point are assigned a value of 0. The sample confidence mask matrix is ​​then backfilled according to the original index.

[0101] The periodic dimension weighted aggregation process is performed on the robust periodic observation matrix to obtain the continuous robust observation matrix. The noise covariance matrix is ​​calculated for the sample intervals with energy envelopes below the threshold. The continuous robust observation matrix is ​​then whitened based on the noise covariance matrix to obtain the preprocessed observation matrix.

[0102] In this embodiment 1, for the continuous robust observation matrix, sample intervals with energy envelopes below a threshold are selected and mean values ​​are removed from each channel to obtain a noise sample matrix; based on the noise sample matrix, the inter-channel covariance is calculated to obtain a noise covariance matrix; eigenvalue decomposition is performed on the noise covariance matrix to obtain an eigenvector matrix and an eigenvalue sequence; the square root reciprocal of the eigenvalues ​​greater than 0 in the eigenvalue sequence is calculated and a diagonal scaling matrix is ​​constructed and multiplied with the eigenvector matrix to obtain a whitening transformation matrix; the continuous robust observation matrix is ​​mean-removed from each channel and left-multiplied by the whitening transformation matrix to obtain a whitened observation matrix as a preprocessed observation matrix.

[0103] In this embodiment, obtaining the purification observation matrix specifically involves:

[0104] The cross-channel energy envelope sequence is calculated based on the preprocessed observation matrix, and autocorrelation analysis is performed to obtain the estimated value of the dominant cardiac cycle and the corresponding delay value range.

[0105] Multiple delay parameters are selected within the delay value range. The covariance matrix of the corresponding delay is calculated for the preprocessed observation matrix. Matrix confidence weights are generated based on the diagonal dominance index of each covariance matrix to obtain a set of multi-delay covariance matrices and a set of matrix confidence weights.

[0106] In this embodiment 1, the diagonal energy is obtained by summing the absolute values ​​of the diagonal elements of the covariance matrix corresponding to each set of delay parameters within the delay value interval; the total energy is obtained by summing the absolute values ​​of all elements of the same covariance matrix; the diagonal dominance index is the ratio of the diagonal energy to the total energy; the total dominance is obtained by summing the diagonal dominance indices corresponding to all delay parameters; the matrix confidence weight is the ratio of the diagonal dominance index to the total dominance.

[0107] The set of multi-delay covariance matrices and the set of matrix confidence weights are used as inputs. The second-order statistical blind source separation algorithm is used to perform weighted joint diagonalization processing to output the separation matrix. The source component sequence is obtained by linear transformation of the preprocessed observation matrix by the separation matrix. At the same time, the equivalent mixing matrix consistent with the separation matrix is ​​calculated.

[0108] In this embodiment 1, an augmented matrix is ​​constructed. The left side of the augmented matrix consists of a separation matrix, and the right side consists of an identity matrix with the same row and column dimensions as the separation matrix. Elementary row operations are performed on the augmented matrix while keeping the row operations synchronously applied to the corresponding columns of the identity matrix until the left side of the augmented matrix is ​​transformed into an identity matrix. The matrix obtained on the right side of the augmented matrix is ​​denoted as the equivalent mixture matrix, thus obtaining the equivalent mixture matrix that satisfies the condition that the product of the separation matrix and the equivalent mixture matrix is ​​an identity matrix.

[0109] For each source component sequence, calculate the set of autocorrelation peak indices within the time delay range;

[0110] In this embodiment 1, each source component sequence is extracted one by one from the source component sequence; the starting delay of the delay value interval is rounded to the minimum delay and the ending delay is rounded to the maximum delay, generating a list containing all integer delays from the minimum to the maximum delay; the mean of the entire sequence is calculated for each source component sequence, and the mean is subtracted from each sample value to obtain the mean-free sequence; for each delay in the integer delay list, the mean-free sequence and the overlapping interval of the mean-free sequence after corresponding delay shift are multiplied point by point and summed, and then divided by the number of sample pairs in the overlapping interval to obtain the corresponding delay autocorrelation value; each delay autocorrelation value is divided by the autocorrelation value with a delay of 0 to obtain the normalized autocorrelation value sequence; the normalized autocorrelation value with the largest value in the normalized autocorrelation value sequence is selected as the autocorrelation peak amplitude, and the delay that produces the maximum value is recorded as the autocorrelation peak delay; the process is repeated for all source component sequences, and the autocorrelation peak amplitude and autocorrelation peak delay are collected to form an autocorrelation peak index set;

[0111] Based on the column vectors of the equivalent mixing matrix, a set of basis vectors for the source space is constructed. The projection energy of each column vector on the basis vector set for the source space and the residual energy on the orthogonal complement space are calculated respectively to obtain a set of spatial projection energy ratio indices.

[0112] Monotonic normalization is performed on the set of autocorrelation peak indices and the set of spatial projection energy ratio indices to generate interference suppression weights;

[0113] The source component sequence is subjected to component-level weighted suppression based on the interference suppression weight, and the weighted source component sequence is reconstructed through an equivalent mixing matrix to output a cleaned observation matrix.

[0114] In this embodiment, the coupling of the separation matrix and the interference suppression weights to generate inversion prior parameters, and the establishment of a weighted forward model based on the inversion prior parameters, specifically involves:

[0115] The source component sequences are sorted according to the interference suppression weights to generate a set of valid source component indices;

[0116] In this embodiment 1, a one-to-one correspondence is established between the interference suppression weights and each component of the source component sequence, forming a pairing table of component numbers and weight values; the pairing table is sorted according to the weight values ​​from largest to smallest to obtain a sorted sequence and a component number mapping sequence; based on the sorted sequence, the ratio of adjacent weight values ​​is calculated and the sorting position number where the ratio reaches the maximum value of the entire sequence is determined. Sort position number The location of the maximum breakpoint in the corresponding weight value sequence; before selecting the mapping sequence. Each component number constitutes a valid source component index set;

[0117] Based on the set of valid source component indices, the corresponding column vectors are extracted from the separation matrix to form a valid separation submatrix. The column vectors of the valid separation submatrix are then scaled to obtain the coupling separation matrix.

[0118] In this embodiment 1, based on the effective source component index set, corresponding weight values ​​are extracted from the interference suppression weights to form an effective weight sequence. The weights of the effective weight sequence are calculated to obtain a normalized benchmark. Column vectors are extracted column by column from the effective separation submatrix based on the effective source component index set to form a column vector sequence. For each column vector in the column vector sequence, the ratio of the corresponding weight value to the normalized benchmark is used as a scaling factor, and the scaling factor is multiplied into each element of the column vector to obtain a scaled column vector. The Euclidean length of the scaled column vector is calculated, and each element of the scaled column vector is divided by the Euclidean length to obtain a normalized column vector. When the Euclidean length is 0, each element of the normalized column vector remains 0. The normalized column vectors are arranged in column order to form a coupling separation matrix.

[0119] The reliability vector of each channel is calculated based on the coupling separation matrix, and the interference suppression weights corresponding to the effective source component index set are monotonically normalized to generate the component confidence sequence.

[0120] In this embodiment 1, the row energy of the coupling separation matrix is ​​calculated according to the row vector corresponding to the channel. The row energy is the sum of the squares of each element of the row vector. The row energy of each channel is divided by the sum of the row energies of all channels to obtain the channel reliability vector. The sum of each element of the channel reliability vector is 1. The effective weight sequence is extracted from the interference suppression weight based on the effective source component index set, and the maximum and minimum values ​​of the effective weight sequence are calculated. Under the condition that the maximum value is greater than the minimum value, the minimum value is subtracted from each weight value of the effective weight sequence to form a difference sequence. Each element of the difference sequence is divided by the difference between the maximum and minimum values ​​to obtain the component confidence sequence. Under the condition that the maximum and minimum values ​​are the same, each element of the component confidence sequence is 1.

[0121] The channel reliability vector and component confidence sequence are combined to generate inversion prior parameters, and the channel weighting factor is generated based on the channel reliability vector, and the source space weighting factor is generated based on the component confidence sequence.

[0122] The spatial coordinates of each channel sensor recorded during multi-channel magnetocardiography acquisition are obtained, and a source space voxel coordinate set is generated within the preset cardiac source space boundary. The mapping coefficient of each voxel to the observation contribution of each channel sensor is calculated based on the spatial coordinates of each channel sensor and the source space voxel coordinate set, forming an unweighted forward model.

[0123] In this embodiment 1, the minimum and maximum values ​​of the spatial coordinates of each channel sensor in the three coordinate axes are calculated based on the channel coordinate table, and the median of the nearest neighbor distance is calculated to generate the sensor bounding box and the median of the nearest neighbor distance. The sensor bounding box is expanded to both sides in the x-axis direction by the median of the nearest neighbor distance, and to both sides in the y-axis direction by the median of the nearest neighbor distance, and then expanded to 6 times the median of the nearest neighbor distance in the negative z-axis direction to form the cardiac source space boundary. Within the cardiac source space boundary, a three-dimensional regular voxel coordinate set is generated by using a grid spacing equal to the median of the nearest neighbor distance divided by 2, and numbered in ascending order of x-coordinate, y-coordinate, and z-coordinate to obtain the source space voxel coordinate set. For each source space voxel coordinate and each sensor space coordinate, the spatial difference vector and Euclidean distance are calculated. The magnetic field vector generated by the unit current dipole in the three orthogonal directions is calculated based on the reciprocal of the cube of the Euclidean distance, and the projection is taken along the sensor sensitive direction to obtain the mapping coefficient. The mapping coefficients corresponding to each source space voxel coordinate are arranged in the order of sensor channel number to form a coefficient vector, and then spliced ​​in the order of voxel number to form an unweighted forward model.

[0124] The channel dimension of the unweighted forward model is weighted according to the channel weighting factor, and the voxel dimension of the unweighted forward model is weighted according to the source space weighting factor. Then, the column vector normalization process is performed on the coefficient vector corresponding to each voxel after weighting to obtain the weighted forward model.

[0125] In this embodiment, the improved FOCUSS algorithm is specifically as follows:

[0126] Under the weighted forward model constraints, the minimum norm solution of the cleaned observation matrix that satisfies the observation fitting constraints is obtained to generate the initial source distribution;

[0127] In this embodiment 1, the purified observation matrix is ​​split into an observation vector sequence according to time samples, and the weighted forward model is read to form a mapping coefficient matrix. For each observation vector, the product of the mapping coefficient matrix and the transpose matrix is ​​calculated to obtain the channel domain Gram matrix, and singular value decomposition is performed on the channel domain Gram matrix. Based on the singular value decomposition result, the Moore-Penrose generalized inverse of the channel domain Gram matrix is ​​constructed and left-multiplied by the observation vector to obtain the channel domain coefficient vector. The transpose matrix of the mapping coefficient matrix is ​​right-multiplied by the channel domain coefficient vector to obtain the source space solution vector, and the residual vector is calculated. The residual vector is obtained by multiplying the mapping coefficient matrix by the source space solution vector and subtracting the observation vector. The source space solution vector satisfies the condition that the source space solution vector has the minimum L2 norm in the solution set where the L2 norm of the residual vector is minimized, and is used as the minimum norm solution. The source space solution vector sequence is obtained by repeating the calculation column by column for the observation vector sequence and then concatenating the sequence column by column to form the initial source distribution.

[0128] Calculate the voxel energy corresponding to each voxel in the source space voxel coordinate set based on the initial source distribution, and establish a three-dimensional adjacency relationship on the source space voxel coordinate set to generate a voxel adjacency table.

[0129] In this embodiment 1, the initial source distribution is arranged into a matrix consisting of voxel row vectors and time sample column vectors according to voxel indices; the voxel energy is obtained by summing the squares of all time sample amplitudes for each voxel row vector and forming a voxel energy vector; the x-axis coordinate sequence, y-axis coordinate sequence, and z-axis coordinate sequence are extracted from the source space voxel coordinate set, and the minimum positive value of the difference between different coordinates is calculated to obtain the x-axis step size, y-axis step size, and z-axis step size; six candidate adjacent coordinates are generated for each voxel coordinate, and the candidate adjacent coordinates differ from the current voxel coordinates by the corresponding step size in one coordinate axis direction and are the same in the other two coordinate axes; the voxel index matching the candidate adjacent coordinates is retrieved in the source space voxel coordinate set, and the matching index is written into the adjacent index list corresponding to the current voxel index to form a voxel adjacency table;

[0130] Select the set of voxels with energy peaks based on voxel energy, and use the set of voxels with energy peaks as the starting point for region growth. Generate a set of connected candidate regions under the constraint of voxel adjacency list according to the expansion order of voxel energy from high to low.

[0131] For each connected candidate domain in the set of connected candidate domains, calculate the domain energy and the number of domain voxels, and calculate the domain-level confidence weights based on the domain energy and the number of domain voxels to obtain the set of domain-level confidence weights;

[0132] In this embodiment 1, for each connected candidate domain, the corresponding voxel energy is extracted from the voxel energy vector according to the voxel index list and summed to obtain the domain energy; for each connected candidate domain, the number of elements in the voxel index list is counted to obtain the domain voxel count; for each connected candidate domain, the quotient of the domain energy and the domain voxel count is calculated to obtain the domain average energy; the domain average energies of all connected candidate domains are summed to obtain the average energy sum; for each connected candidate domain, the domain average energy is divided by the average energy sum to obtain the domain-level confidence weight; the domain-level confidence weights of all connected candidate domains are arranged in the order of connected candidate domains to form a domain-level confidence weight set;

[0133] A connected-domain focused sparse weight vector is generated based on the domain-level confidence weight set, and the voxel dimension of the weighted forward model is scaled to obtain the connected-domain weighted forward model.

[0134] In this embodiment 1, the connected candidate domain set, the domain-level confidence weight set, the source space voxel coordinate set, and the weighted forward model are read to establish a correspondence between the column index of the weighted forward model and the voxel index of the source space voxel coordinate set. The connected domain focused sparse weight vector is initialized, with its length equal to the number of voxels in the source space voxel coordinate set and each element having a value of 0. The connected candidate domain set is traversed, and the elements of the domain-level confidence weight set are read as domain weights according to the traversal sequence. The domain weights are written to the corresponding voxel index positions of the connected domain focused sparse weight vector according to the connected candidate domain voxel index list. When a voxel index position is written repeatedly, the maximum value of the domain weight is retained. All voxel column vectors of the weighted forward model are traversed, and the connected domain focused sparse weight vector is multiplied by the index element of the corresponding voxel column vector to obtain a scaled column vector, which is then used to replace the original voxel column vector to obtain the connected domain weighted forward model.

[0135] Under the constraints of the connected component weighted forward model, the weighted least squares subproblem of cleaning the observation matrix is ​​solved to obtain the intermediate solution vector. The intermediate solution vector is then merged and updated with the connected component focused sparse weight vector to generate the current source distribution.

[0136] The residual matrix is ​​calculated based on the weighted forward model of connected components and the current source distribution. The ratio of the number of elements in the intersection to the number of elements in the union of the candidate connected component sets in two adjacent rounds is also calculated. The energy change of the residual matrix in two adjacent rounds is calculated to generate the convergence determination quantity.

[0137] In this embodiment 1, the predicted observation matrix is ​​obtained by multiplying the weighted forward model of the connected domains with the current source distribution; the residual matrix is ​​obtained by subtracting the purified observation matrix from the predicted observation matrix; the residual energy of the current round is obtained by summing the squares of the amplitudes of all elements of the residual matrix; the residual energy of the previous round is read, the absolute value of the difference between the residual energy of the current round and the residual energy of the previous round is taken, and the ratio is calculated with the residual energy of the previous round to obtain the change in residual energy; the connected candidate domain set of the current round and the connected candidate domain set of the previous round are read, and all voxel indices in each connected candidate domain set are deduplicated and merged to obtain the voxel set of the current round and the voxel set of the previous round; the number of intersection voxels and the number of union voxels are calculated for the voxel set of the current round and the voxel set of the previous round, and the intersection-union ratio is calculated by ratio of the number of intersection voxels and the number of union voxels; the domain change is obtained by subtracting the intersection-union ratio from 1, and the domain change is summed with the residual energy change to obtain the convergence determination value;

[0138] When the convergence determination quantity meets the preset termination threshold, the convergence source distribution is output. When the convergence determination quantity does not meet the preset termination threshold, the current source distribution is used as the initial source distribution for the next round and the iteration count is updated. When the maximum number of iterations is reached, the convergence source distribution is output.

[0139] In this embodiment 1, a preset termination threshold is set. Specifically, the convergence determination quantity is obtained by summing the domain change quantity and the residual energy change quantity, with the domain change quantity ranging from 0 to 1 and the residual energy change quantity ranging from 0 to 1. The termination threshold is set to 0.05. When the convergence determination quantity is less than or equal to 0.05, it is determined that the stability of the intersection-union ratio of the connected candidate domain set and the energy stability of the residual matrix satisfy the termination condition, and the converged source distribution is output. When the convergence determination quantity is greater than 0.05, the current source distribution is iteratively updated and the iteration count is updated.

[0140] In this embodiment, the localization process performed on the convergence source distribution specifically includes:

[0141] Perform voxel energy normalization on the convergence source distribution, calculate the voxel energy statistics of the convergence source distribution, and generate energy quantile threshold parameters and normalized energy voxel maps.

[0142] In this embodiment 1, the source strength value of the convergent source distribution corresponding to each voxel in the source space voxel coordinate set is read, and the source strength value is squared to generate a voxel energy vector and a voxel energy voxel map is generated simultaneously; the voxel energy minimum, voxel energy maximum, and voxel energy median value are calculated to form voxel energy statistics; the voxel energy vector is sorted in ascending order, and the product of the voxel energy vector length and 0.90 is rounded up to obtain the quantile number, and the voxel energy corresponding to the quantile number is taken as the energy quantile threshold parameter; the voxel energy voxel map is linearly scaled using the voxel energy minimum and voxel energy maximum as boundaries to obtain a normalized energy voxel map, and the normalized energy voxel map is assigned a value of 0 when the voxel energy maximum equals the voxel energy minimum;

[0143] Voxels whose voxel energies are not lower than the energy quantile threshold parameter in the normalized energy voxel graph are marked as high-energy voxels. A three-dimensional six-adjacency relationship is established on the source space voxel coordinate set, and a set of peak connected regions is generated based on the high-energy voxels.

[0144] In this embodiment 1, voxels marked as 1 are selected from the high-energy voxel marking results to form a set of voxels to be searched, and an unvisited flag vector corresponding to the set of voxels to be searched is generated; unvisited voxels are extracted from the set of voxels to be searched as seed voxels, and a breadth-first traversal is performed based on the three-dimensional six-adjacency relationship. All connected high-energy voxels obtained by the traversal are written into the connected voxel set and the unvisited flag vector is updated; the voxel energy of the connected voxel set is extracted from the normalized energy voxel graph and the maximum voxel energy is calculated. The voxel corresponding to the maximum voxel energy is recorded as the peak voxel; the connected voxel set and the peak voxel form a peak connected domain, which is assigned a domain number and added to the peak connected domain set, until there are no voxels marked as 1 in the unvisited flag vector;

[0145] Calculate the domain confidence weight set based on the peak connected component set;

[0146] In this embodiment 1, for each peak connected domain, the energy values ​​of all voxels in the connected voxel set are extracted from the normalized energy voxel graph, the domain energy sum is calculated, and the number of domain voxels is calculated. Under the constraint of three-dimensional six adjacency relations, the number of boundary voxels in the connected voxel set is counted. Boundary voxels satisfy that at least one of their six adjacent voxels does not belong to the connected voxel set. The energy concentration is obtained by dividing the voxel energy corresponding to the peak voxel by the domain energy sum, and the energy concentration is set to 0 when the domain energy sum equals 0. The shape compactness is obtained by subtracting the number of boundary voxels from the number of domain voxels and dividing it by the number of domain voxels, and the shape compactness is set to 0 when the number of domain voxels equals 0. The minimum and maximum values ​​of the energy concentration and shape compactness of all peak connected domains are calculated respectively, and linear normalization is performed to obtain the normalized concentration and normalized compactness. When the maximum value equals the minimum value, the corresponding normalization result is set to 1. The normalized concentration and normalized compactness are multiplied to obtain the domain confidence weight, and the domain confidence weight set is output according to the peak connected domain number.

[0147] Based on the set of domain confidence weights, peak connected regions with domain confidence weights greater than a preset threshold are selected. The selected peak connected regions are merged to obtain a set of candidate lesion voxels. A layer of adjacent voxels is then extended outside the set of candidate lesion voxels to generate a boundary search voxel band.

[0148] The spatial energy gradient vector field of the convergent source distribution is calculated on the boundary search voxel band, and non-maximum suppression is performed to obtain the candidate boundary voxel set.

[0149] In this embodiment 1, the energy values ​​of the forward and reverse adjacent voxels in the coordinate axes 1, 2, and 3 directions are extracted voxels one by one within the boundary search voxel band. The difference between the energy values ​​of the forward and reverse adjacent voxels is divided by the corresponding coordinate difference to obtain the gradient components in three directions, which are then used to form a gradient vector. At the same time, the magnitude of the gradient vector is calculated to obtain the gradient magnitude voxel map. For the gradient vector of each voxel, the absolute values ​​of the gradient components in three directions are calculated, and the direction corresponding to the largest absolute value is selected as the gradient principal direction. The forward and reverse adjacent voxels are determined based on the gradient principal direction. The gradient magnitude of the current voxel is compared with the gradient magnitudes of the forward and reverse adjacent voxels. When the gradient magnitude of the current voxel is greater than both the gradient magnitudes of the forward and reverse adjacent voxels, the current voxel is marked. All voxels marked as 1 are output to form a candidate boundary voxel set.

[0150] The candidate boundary voxel set is subjected to closed connectivity verification, boundary segments that do not form closed toruses are deleted, and the gap voxels are subjected to connectivity completion processing to obtain the lesion voxel boundary set.

[0151] In this embodiment, the three-dimensional voxel aggregation process performed on the lesion voxel boundary set specifically includes:

[0152] Map the set of voxel boundaries of the lesion to the set of voxel coordinates of the source space to generate a boundary occupation marker voxel map and a set of boundary voxels.

[0153] Based on the three-dimensional six-adjacency relationship of the source space voxel coordinate set, the boundary voxel set is subjected to connectivity aggregation processing, and the parity counting and internal / external discrimination are performed along three orthogonal directions to obtain the lesion entity voxel set;

[0154] In this embodiment 1, three mutually orthogonal coordinate axes are defined based on the source space voxel coordinate set. For each voxel to be identified within the boundary-occupying marker voxel map, starting from the center of the voxel to be identified, the ray path is moved along the positive direction of the first coordinate axis to the outside of the source space voxel coordinate set, and the number of times the ray path crosses the boundary voxel set is recorded. Voxels to be identified with an odd number of crossings are marked as inner markers, and voxels to be identified with an even number of crossings are marked as outer markers, thus obtaining the first direction identification result. The second and third direction identification results are obtained using the same voxel-by-voxel stepping and counting method. When the number of inner markers in the three direction identification results is not less than 2, the voxel to be identified is determined to be located inside the lesion entity voxel set and added to the lesion entity voxel set. When the number of inner markers in the three direction identification results is less than 2, the voxel to be identified is determined to be located outside the lesion entity voxel set and removed from the lesion entity voxel set.

[0155] Calculate the discrete shortest distance from the voxel to the boundary occupied marker voxel map based on the voxel set of lesion entities, generate a discrete distance voxel map, and extract the candidate voxel set of holes from the discrete distance voxel map.

[0156] In this embodiment 1, voxels that do not belong to the lesion entity voxel set are enumerated based on the source space voxel coordinate set and recorded as the background voxel set; within the background voxel set, background voxels located at the outer boundary of the source space voxel coordinate set are selected as seed voxels and connected expansion is performed based on the three-dimensional six-adjacency relationship to obtain the outer background voxel set; the outer background voxel set is removed from the background voxel set to obtain the inner background voxel set; within the inner background voxel set, voxels with a discrete distance voxel map distance value of 0 are removed and voxels with a distance value greater than 0 are retained to obtain the hole candidate voxel set;

[0157] Based on the discrete distance voxel map, structural element scale parameters are generated at the corresponding positions of the candidate voxel set of pores, and three-dimensional morphological closing operations are performed on the lesion entity voxel set to obtain the pore repair voxel set.

[0158] The sign distance value of the voxel corner points is calculated based on the voxel set of hole repair, and the intersection point is determined at the voxel boundary to obtain the initial three-dimensional lesion surface mesh;

[0159] The non-manifold edges and boundary loops of the initial 3D lesion surface mesh are detected, and boundary loop filling is performed. Local Laplacian smoothing is then applied to the filled area to output the 3D lesion surface mesh.

[0160] In this embodiment, the consistency verification process between the three-dimensional lesion surface mesh and the convergence source distribution is specifically as follows:

[0161] The coordinates of the grid vertices of the three-dimensional lesion surface mesh are transformed to the source space voxel coordinate set coordinate system. For each grid vertex, the voxel unit containing the grid vertex is determined, and a neighborhood voxel set and a vertex sampling index table are generated.

[0162] The voxel energy of the neighborhood voxel set in the convergence source distribution is read from the vertex sampling index table, and a weighted summation is performed according to the distance attenuation weight from the voxel center to the grid vertex to generate an energy projection value sequence.

[0163] Calculate the projection value statistics and projection value quantile threshold parameters based on the energy projection value sequence. Mark the grid vertices whose projection value statistics are lower than the projection value quantile threshold parameters as low-consistency vertices to obtain a set of low-consistency vertices.

[0164] In this embodiment 1, the energy projection value sequence is mapped to the mesh vertices one by one according to the vertex sampling index table; the set of adjacent vertices sharing edges with the target mesh vertex is extracted from the triangular patch set of the three-dimensional lesion surface mesh; the projection values ​​corresponding to the target mesh vertex and the adjacent vertices are read to form a local projection value set; the local projection value set is sorted in ascending order and the median is taken as the projection value statistic of the target mesh vertex; the projection value statistics are formed into a statistical sequence and sorted in ascending order; the product of 0.25 and the sequence length is calculated to obtain the index m, and the index m is taken as the smallest integer not less than the product; the m-th value is taken as the projection value quantile threshold parameter; the mesh vertices whose projection value statistics are less than the projection value quantile threshold parameter are marked as low consistency vertices;

[0165] In the 3D lesion surface mesh, delete the triangles adjacent to the set of low-consistency vertices while maintaining the connectivity of the remaining triangles to generate a trimmed lesion surface mesh.

[0166] Based on the mesh of the lesion surface, a mesh-occupied voxel map is constructed, and voxels whose voxel energy is not lower than the energy quantile threshold parameter in the convergence source distribution are marked as high-energy voxel sets. The proportion of voxels falling outside the mesh-occupied voxel map of the high-energy voxel sets is calculated to generate an external leakage consistency index.

[0167] In this embodiment 1, the vertex coordinates of the triangular pieces of the lesion surface mesh are read and transformed to the source space voxel coordinate set coordinate system; for each triangular piece, the minimum voxel index and the maximum voxel index corresponding to the vertex coordinates are calculated to form a voxel bounding box; all voxel units within the voxel bounding box are traversed, and the overlap judgment of the voxel unit cube and the triangular piece is performed based on the separating axis theorem to obtain the overlap mark; the voxel units with the overlap mark corresponding to the true voxel are written into the mesh-occupied voxel map and assigned a value of 1, while the voxel units not written into the voxel are kept with a value of 0; after all triangular pieces are processed, the mesh-occupied voxel map is output.

[0168] When the leakage consistency index is higher than the preset leakage threshold, where the preset leakage threshold is set to 0.10, the triangles adjacent to the leaked voxels outside the grid are deleted, and boundary ring filling is performed to obtain the three-dimensional lesion model.

[0169] In this embodiment, a lesion model construction system based on magnetocardiography includes:

[0170] The multi-channel acquisition and coordinate recording module acquires raw observation data of multi-channel magnetocardiogram and records sensor spatial coordinates, sampling rate parameters and channel number information, and outputs raw observation matrix and channel coordinate table;

[0171] The preprocessing module performs spatial basis function fitting to subtract far-field interference, band-limited coherence weighting, heart rate cycle alignment, and noise covariance whitening on the original observation matrix, and outputs the preprocessed observation matrix.

[0172] The weighted joint diagonalization purification module constructs a set of multi-delay covariance matrices and a set of matrix confidence weights, performs weighted joint diagonalization processing, and outputs a purified observation matrix.

[0173] The inversion prior coupling and weighted forward modeling module performs coupling on the separation matrix and interference suppression weights to generate inversion prior parameters, generates channel weighting factors and source space weighting factors, and establishes a weighted forward model;

[0174] The connected domain focusing FOCUSS sparse reconstruction module generates a set of connected candidate domains based on voxel energy and three-dimensional adjacency relations, calculates domain-level confidence weights, constructs a connected domain weighted forward model, and iteratively solves for the output convergence source distribution.

[0175] The focalization and boundary extraction module generates a set of peak connected regions and a set of candidate lesion voxels based on the energy quantile threshold parameter, calculates the spatial energy gradient vector field and performs non-maximum suppression and closed connectivity verification, and outputs a set of lesion voxel boundaries.

[0176] The mesh reconstruction and consistency verification module performs three-dimensional voxel aggregation, morphological closing operation and surface triangulation reconstruction on the lesion voxel boundary set to obtain a three-dimensional lesion surface mesh, calculates the energy projection value sequence and leakage consistency index and trims and repairs it, and outputs a three-dimensional lesion model.

[0177] Example 2:

[0178] To verify the feasibility of this invention in practice, it was applied to a real-world scenario of multi-channel magnetocardiography (MCG) lesion localization and 3D lesion model construction. Five sets of MCG data were selected and a corresponding reference 3D lesion model was obtained to test the localization stability under conditions of strong noise interference and significant individual differences. In this scenario, existing methods are prone to instability in the source component sequence when preprocessing the statistical structure fluctuations of the observation matrix, which in turn leads to discrete sparse point pseudo-lesions in sparse reconstruction and causes surface mesh leakage. This invention obtains the convergent source distribution by coupling interference suppression weights with inversion prior parameters and using an improved FOCUSS algorithm with connected component focusing.

[0179] During implementation, raw magnetocardiogram (MCC) data from 64 channels was collected, and sensor spatial coordinates, sampling rate parameters, and channel numbers were recorded. The sampling rate was 1000 Hz, and the acquisition time was 20 seconds, forming a raw observation matrix and a channel coordinate table. Spatial basis function fitting was performed on the raw observation matrix to subtract far-field interference, resulting in a far-field-de-field observation matrix. Band-limited coherence indices were calculated in the 5 Hz to 45 Hz range, and channel reliability weight vectors were generated and applied to obtain a reliability-weighted observation matrix. Based on the reliability-weighted observation matrix, cardiac marker sequences were calculated, and periodic slicing and periodic resampling were performed to align them, resulting in an aligned periodic observation matrix. Based on the aligned periodic observation matrix, a sample confidence mask matrix was calculated and applied to obtain a robust periodic observation matrix. Periodic dimension weighted aggregation was performed to obtain a continuous robust observation matrix. Sample intervals with energy envelopes below the 0.20 quantile were selected, and the noise covariance matrix was calculated and whitened to obtain a preprocessed observation matrix.

[0180] In the construction of the purified observation matrix, a set of multi-delay covariance matrices is calculated based on the preprocessed observation matrix, and a set of matrix confidence weights is generated. Weighted joint diagonalization is performed to output a separation matrix. The source component sequences and equivalent mixing matrices are obtained by linearly transforming the preprocessed observation matrix using the separation matrix. The set of autocorrelation peak indices and spatial projection energy ratio indices are calculated and monotonically normalized to generate interference suppression weights. Component-level weighted suppression is applied to the source component sequences based on these interference suppression weights, and the purified observation matrix is ​​reconstructed using the equivalent mixing matrix. The source component sequences are sorted according to the interference suppression weights to generate a set of effective source component indices. The corresponding column vectors of the separation matrix are extracted to form effective separation sub-matrices, and column vector scaling is performed to obtain the coupling separation matrix. The channel reliability vector is calculated, and a component confidence sequence is generated. These are combined to obtain the inversion prior parameters and generate channel weighting factors and source spatial weighting factors. Within the preset cardiac source space boundary, a set of voxel coordinates with a voxel side length of 3 mm is generated, with 21,000 voxels. The mapping coefficients are calculated to form an unweighted forward model, and then weighted processing of channel dimension and voxel dimension and column vector normalization processing are applied to obtain a weighted forward model.

[0181] In the improved FOCUSS algorithm, the minimum norm solution of the cleaned observation matrix is ​​first solved to generate the initial source distribution. Voxel energies are calculated, and a three-dimensional adjacency table is generated by establishing adjacency relationships. A set of voxels with peak energy is selected to generate a set of connected candidate domains, and a set of domain-level confidence weights is calculated. A connected domain focusing sparse weight vector is generated, and the weighted forward model is scaled to obtain a connected domain weighted forward model. The weighted least squares subproblem is solved to obtain intermediate solution vectors, which are then merged and updated to obtain the current source distribution. A preset termination threshold is set with an intersection-union ratio (IU) of at least 0.95 and a residual matrix energy change of at least 0.01. The maximum number of iterations is 30, and the convergence iterations for the 5 sets of samples are 18, 16, 19, 17, and 15, outputting the converged source distribution.

[0182] The convergence source distribution is localized and an energy quantile threshold parameter of 0.85 is generated. After obtaining the lesion voxel boundary set, three-dimensional voxel aggregation is performed to obtain the pore repair voxel set. The structural element scale parameter is limited to 2 to 4 voxels, and surface triangulation reconstruction is completed to output the three-dimensional lesion surface mesh. When performing consistency verification processing on the three-dimensional lesion surface mesh and the convergence source distribution, the neighborhood voxel set is taken as the voxel unit and its 26 neighboring voxels. The distance attenuation weight adopts the inverse square form, the projection value quantile threshold parameter is set to 0.25 quantile, and the leakage threshold is set to 0.12 to obtain the three-dimensional lesion model. The number of triangular pieces in the three-dimensional lesion surface mesh is 8200 to 11500, the number of local Laplacian smoothing iterations is 12, and the calculation time for a single sample is 8.4 seconds to 11.2 seconds.

[0183] After adopting the scheme of this invention, the voxel overlap coefficient reached 0.82 to 0.88, the lesion center localization error was 2.1 mm to 4.3 mm, and the leakage consistency index was no higher than 0.10. In contrast, the control scheme, without using the matrix confidence weight set and connected component focusing sparse weight vector, had a voxel overlap coefficient of 0.71 to 0.79, a lesion center localization error of 5.8 mm to 9.6 mm, and a leakage consistency index as high as 0.19. The table below shows a comparison between the reference measured values ​​and the predicted values ​​of this invention for five groups of samples.

[0184] Table 1 Comparison of Measured and Predicted Values ​​of Geometric Indicators for Lesion Models

[0185] Sample number Measured distance (in millimeters) from the center of the lesion to the reference point Predicted distance in millimeters from lesion center to reference point The maximum diameter of the lesion was measured in millimeters. Predicted maximum diameter of lesion (in millimeters) Measured lesion volume (cubic millimeters) Lesion volume prediction (cubic millimeters) Actual measured surface area of ​​lesion (square millimeters) Predicted lesion surface area in square millimeters Actual measurement of leakage consistency index External leakage consistency index prediction MCG-001 48.6 50.8 26.4 27.2 6120 6460 2560 2705 0.06 0.07 MCG-002 54.1 57.2 30.8 31.9 7420 7810 2985 3150 0.08 0.09 MCG-003 41.9 44.6 21.6 22.2 3980 4140 1960 2050 0.05 0.06 MCG-004 58.7 62.5 32.1 33.5 8120 8640 3220 3408 0.09 0.10 MCG-005 46.3 48.7 24.9 25.8 5340 5570 2340 2452 0.07 0.08

[0186] As shown in Table 1, the predicted and measured distances from the lesion center to the reference point remained within the range of 2.2 mm to 3.8 mm, and the maximum diameter difference remained within the range of 0.6 mm to 1.4 mm, indicating that the convergence source distribution can stably converge to a single peak connected domain after focalization processing. The relative error of the predicted lesion volume was no higher than 6.5%, and the relative error of the predicted lesion surface area was no higher than 5.9%, indicating that the voids in the lesion voxel boundary set were sufficiently repaired after 3D morphological closing operations, and the 3D lesion surface mesh obtained by surface triangulation reconstruction was consistent with the geometric scale of the reference model. The predicted values ​​of the leakage consistency index were all below 0.10 and close to the measured values. Combined with the removal of low consistency vertices and leakage threshold constraints in the consistency verification process, it can be seen that the proportion of high-energy voxels on the outside of the voxel map occupied by the mesh was effectively suppressed, and the leakage risk of the 3D lesion model was reduced. The control scheme exhibits scattered points in multiple connected components under the same data, leading to an increase in the leakage consistency index. This aligns with the mechanism of this invention, which uses matrix confidence weight sets to stably and jointly diagonalize the output and shrinks the support set by focusing sparse weight vectors on connected components. Furthermore, the improved FOCUSS algorithm reaches the preset termination threshold within 15 to 19 iterations in 5 sets of samples, with residual matrix energy changes all below 0.01, indicating that the convergence criterion exhibits stable termination behavior and reduces redundant computation overhead.

[0187] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for constructing a lesion model based on magnetocardiography, characterized in that, include: Raw observation data from multi-channel magnetocardiography were collected and preprocessed to obtain a preprocessed observation matrix; The set of multi-delay covariance matrices is calculated based on the preprocessed observation matrix, and joint diagonalization is performed by the second-order statistical blind source separation algorithm to obtain the separation matrix and source component sequence. Interference suppression weights are generated based on the autocorrelation peak and spatial projection energy ratio of the source component sequence to obtain the purified observation matrix. The separation matrix and the interference suppression weights are coupled to generate inversion prior parameters, and a weighted forward model is established based on the inversion prior parameters. Under the constraints of the weighted forward model, the improved FOCUSS algorithm is used to perform sparse reconstruction processing on the cleaned observation matrix to obtain the convergence source distribution; The convergence source distribution is localized, and a set of candidate lesion voxels is generated based on the peak connected domain of the convergence source distribution. The boundary of the candidate lesion voxel set is determined based on the energy gradient of the convergence source distribution, and the lesion voxel boundary set is obtained. Three-dimensional voxel aggregation processing is performed on the voxel boundary set of the lesion, morphological closing operation is used to eliminate holes and surface triangulation reconstruction is performed to obtain the three-dimensional lesion surface mesh; Perform consistency verification on the three-dimensional lesion surface mesh and the convergence source distribution, calculate the energy projection value of the source distribution corresponding to the mesh vertex and remove mesh vertices below the threshold, and output the three-dimensional lesion model.

2. The method for constructing a lesion model based on magnetocardiography according to claim 1, characterized in that, The preprocessing process, resulting in a preprocessed observation matrix, is as follows: Collect raw observation data of multi-channel magnetocardiography and simultaneously record the spatial coordinates of each channel sensor, sampling rate parameters and channel number information to obtain the raw observation matrix and channel coordinate table; Based on the channel coordinate table, a spatial basis function is constructed. For each sampling time of the original observation matrix, the spatial basis function is fitted on the entire channel. The fitting result is then subtracted from the original observation matrix to obtain the far-field observation matrix. In the far-field observation matrix, the band-limited coherence index between each channel and its neighboring channel set is calculated. Based on the band-limited coherence index, a channel reliability weight vector is generated and applied to the far-field observation matrix to obtain the reliability-weighted observation matrix. The heartbeat marker sequence is calculated based on the reliability-weighted observation matrix, and the reliability-weighted observation matrix is ​​aligned by periodic slicing and periodic resampling to obtain the aligned periodic observation matrix. The sample confidence mask matrix is ​​calculated based on the aligned periodic observation matrix and applied to the aligned periodic observation matrix to obtain the robust periodic observation matrix; The periodic dimension weighted aggregation process is performed on the robust periodic observation matrix to obtain the continuous robust observation matrix. The noise covariance matrix is ​​calculated for the sample intervals with energy envelopes below the threshold. The continuous robust observation matrix is ​​then whitened based on the noise covariance matrix to obtain the preprocessed observation matrix.

3. The method for constructing a lesion model based on magnetocardiography according to claim 1, characterized in that, The process of obtaining the purified observation matrix is ​​as follows: The cross-channel energy envelope sequence is calculated based on the preprocessed observation matrix, and autocorrelation analysis is performed to obtain the estimated value of the dominant cardiac cycle and the corresponding delay value range. Multiple delay parameters are selected within the delay value range. The covariance matrix of the corresponding delay is calculated for the preprocessed observation matrix. Matrix confidence weights are generated based on the diagonal dominance index of each covariance matrix to obtain a set of multi-delay covariance matrices and a set of matrix confidence weights. The set of multi-delay covariance matrices and the set of matrix confidence weights are used as inputs. The second-order statistical blind source separation algorithm is used to perform weighted joint diagonalization processing to output the separation matrix. The source component sequence is obtained by linear transformation of the preprocessed observation matrix by the separation matrix. At the same time, the equivalent mixing matrix consistent with the separation matrix is ​​calculated. For each source component sequence, calculate the set of autocorrelation peak indices within the time delay range; Based on the column vectors of the equivalent mixing matrix, a set of basis vectors for the source space is constructed. The projection energy of each column vector on the basis vector set for the source space and the residual energy on the orthogonal complement space are calculated respectively to obtain a set of spatial projection energy ratio indices. Monotonic normalization is performed on the set of autocorrelation peak indices and the set of spatial projection energy ratio indices to generate interference suppression weights; The source component sequence is subjected to component-level weighted suppression based on the interference suppression weight, and the weighted source component sequence is reconstructed through an equivalent mixing matrix to output a cleaned observation matrix.

4. The method for constructing a lesion model based on magnetocardiography according to claim 1, characterized in that, The process involves coupling the separation matrix with the interference suppression weights to generate inversion prior parameters, and then establishing a weighted forward model based on these parameters. Specifically: The source component sequences are sorted according to the interference suppression weights to generate a set of valid source component indices; Based on the set of valid source component indices, the corresponding column vectors are extracted from the separation matrix to form a valid separation submatrix. The column vectors of the valid separation submatrix are then scaled to obtain the coupling separation matrix. The reliability vector of each channel is calculated based on the coupling separation matrix, and the interference suppression weights corresponding to the effective source component index set are monotonically normalized to generate the component confidence sequence. The channel reliability vector and component confidence sequence are combined to generate inversion prior parameters, and the channel weighting factor is generated based on the channel reliability vector, and the source space weighting factor is generated based on the component confidence sequence. The spatial coordinates of each channel sensor recorded during multi-channel magnetocardiography acquisition are obtained, and a source space voxel coordinate set is generated within the preset cardiac source space boundary. The mapping coefficient of each voxel to the observation contribution of each channel sensor is calculated based on the spatial coordinates of each channel sensor and the source space voxel coordinate set, forming an unweighted forward model. The channel dimension of the unweighted forward model is weighted according to the channel weighting factor, and the voxel dimension of the unweighted forward model is weighted according to the source space weighting factor. Then, the column vector normalization process is performed on the coefficient vector corresponding to each voxel after weighting to obtain the weighted forward model.

5. The method for constructing a lesion model based on magnetocardiography according to claim 1, characterized in that, The improved FOCUSS algorithm is specifically as follows: Under the weighted forward model constraints, the minimum norm solution of the cleaned observation matrix that satisfies the observation fitting constraints is obtained to generate the initial source distribution; Calculate the voxel energy corresponding to each voxel in the source space voxel coordinate set based on the initial source distribution, and establish a three-dimensional adjacency relationship on the source space voxel coordinate set to generate a voxel adjacency table. Select the set of voxels with energy peaks based on voxel energy, and use the set of voxels with energy peaks as the starting point for region growth. Generate a set of connected candidate regions under the constraint of voxel adjacency list according to the expansion order of voxel energy from high to low. For each connected candidate domain in the set of connected candidate domains, calculate the domain energy and the number of domain voxels, and calculate the domain-level confidence weights based on the domain energy and the number of domain voxels to obtain the set of domain-level confidence weights; A connected-domain focused sparse weight vector is generated based on the domain-level confidence weight set, and the voxel dimension of the weighted forward model is scaled to obtain the connected-domain weighted forward model. Under the constraints of the connected component weighted forward model, the weighted least squares subproblem of cleaning the observation matrix is ​​solved to obtain the intermediate solution vector. The intermediate solution vector is then merged and updated with the connected component focused sparse weight vector to generate the current source distribution. The residual matrix is ​​calculated based on the weighted forward model of connected components and the current source distribution. The ratio of the number of elements in the intersection to the number of elements in the union of the candidate connected component sets in two adjacent rounds is also calculated. The energy change of the residual matrix in two adjacent rounds is calculated to generate the convergence determination quantity. When the convergence determination quantity meets the preset termination threshold, the convergence source distribution is output. When the convergence determination quantity does not meet the preset termination threshold, the current source distribution is used as the initial source distribution for the next round and the iteration count is updated. When the maximum number of iterations is reached, the convergence source distribution is output.

6. The method for constructing a lesion model based on magnetocardiography according to claim 1, characterized in that, The process of performing focalization on the convergence source distribution specifically involves: Perform voxel energy normalization on the convergence source distribution, calculate the voxel energy statistics of the convergence source distribution, and generate energy quantile threshold parameters and normalized energy voxel maps. Voxels whose voxel energies are not lower than the energy quantile threshold parameter in the normalized energy voxel graph are marked as high-energy voxels. A three-dimensional six-adjacency relationship is established on the source space voxel coordinate set, and a set of peak connected regions is generated based on the high-energy voxels. Calculate the domain confidence weight set based on the peak connected component set; Based on the set of domain confidence weights, peak connected regions with domain confidence weights greater than a preset threshold are selected. The selected peak connected regions are merged to obtain a set of candidate lesion voxels. A layer of adjacent voxels is then extended outside the set of candidate lesion voxels to generate a boundary search voxel band. The spatial energy gradient vector field of the convergent source distribution is calculated on the boundary search voxel band, and non-maximum suppression is performed to obtain the candidate boundary voxel set. The candidate boundary voxel set is subjected to closed connectivity verification, boundary segments that do not form closed toruses are deleted, and the gap voxels are subjected to connectivity completion processing to obtain the lesion voxel boundary set.

7. The method for constructing a lesion model based on magnetocardiography according to claim 1, characterized in that, The specific steps of performing three-dimensional voxel aggregation on the lesion voxel boundary set are as follows: Map the set of voxel boundaries of the lesion to the set of voxel coordinates of the source space to generate a boundary occupation marker voxel map and a set of boundary voxels. Based on the three-dimensional six-adjacency relationship of the source space voxel coordinate set, the boundary voxel set is subjected to connectivity aggregation processing, and the parity counting and internal / external discrimination are performed along three orthogonal directions to obtain the lesion entity voxel set; Calculate the discrete shortest distance from the voxel to the boundary occupied marker voxel map based on the voxel set of lesion entities, generate a discrete distance voxel map, and extract the candidate voxel set of holes from the discrete distance voxel map. Based on the discrete distance voxel map, structural element scale parameters are generated at the corresponding positions of the candidate voxel set of pores, and three-dimensional morphological closing operations are performed on the lesion entity voxel set to obtain the pore repair voxel set. The sign distance value of the voxel corner points is calculated based on the voxel set of hole repair, and the intersection point is determined at the voxel boundary to obtain the initial three-dimensional lesion surface mesh; The non-manifold edges and boundary loops of the initial 3D lesion surface mesh are detected, and boundary loop filling is performed. Local Laplacian smoothing is then applied to the filled area to output the 3D lesion surface mesh.

8. The method for constructing a lesion model based on magnetocardiography according to claim 1, characterized in that, The consistency verification process for the three-dimensional lesion surface mesh and the convergence source distribution is specifically as follows: The coordinates of the grid vertices of the three-dimensional lesion surface mesh are transformed to the source space voxel coordinate set coordinate system. For each grid vertex, the voxel unit containing the grid vertex is determined, and a neighborhood voxel set and a vertex sampling index table are generated. The voxel energy of the neighborhood voxel set in the convergence source distribution is read from the vertex sampling index table, and a weighted summation is performed according to the distance attenuation weight from the voxel center to the grid vertex to generate an energy projection value sequence. Calculate the projection value statistics and projection value quantile threshold parameters based on the energy projection value sequence. Mark the grid vertices whose projection value statistics are lower than the projection value quantile threshold parameters as low-consistency vertices to obtain a set of low-consistency vertices. In the 3D lesion surface mesh, delete the triangles adjacent to the set of low-consistency vertices while maintaining the connectivity of the remaining triangles to generate a trimmed lesion surface mesh. Based on the mesh of the lesion surface, a mesh-occupied voxel map is constructed, and voxels whose voxel energy is not lower than the energy quantile threshold parameter in the convergence source distribution are marked as high-energy voxel sets. The proportion of voxels falling outside the mesh-occupied voxel map of the high-energy voxel sets is calculated to generate an external leakage consistency index. When the leakage consistency index is higher than the preset leakage threshold, delete the triangular pieces adjacent to the leaked voxels outside the grid on the voxel map, and perform boundary ring filling to obtain a three-dimensional lesion model.

9. A lesion model construction system based on magnetocardiography, comprising the lesion model construction method based on magnetocardiography as described in any one of claims 1 to 8, characterized in that, include: The multi-channel acquisition and coordinate recording module acquires raw observation data of multi-channel magnetocardiogram and records sensor spatial coordinates, sampling rate parameters and channel number information, and outputs raw observation matrix and channel coordinate table; The preprocessing module performs spatial basis function fitting to subtract far-field interference, band-limited coherence weighting, heart rate cycle alignment, and noise covariance whitening on the original observation matrix, and outputs the preprocessed observation matrix. The weighted joint diagonalization purification module constructs a set of multi-delay covariance matrices and a set of matrix confidence weights, performs weighted joint diagonalization processing, and outputs a purified observation matrix. The inversion prior coupling and weighted forward modeling module performs coupling on the separation matrix and interference suppression weights to generate inversion prior parameters, generates channel weighting factors and source space weighting factors, and establishes a weighted forward model; The connected domain focusing FOCUSS sparse reconstruction module generates a set of connected candidate domains based on voxel energy and three-dimensional adjacency relations, calculates domain-level confidence weights, constructs a connected domain weighted forward model, and iteratively solves for the output convergence source distribution. The focalization and boundary extraction module generates a set of peak connected regions and a set of candidate lesion voxels based on the energy quantile threshold parameter, calculates the spatial energy gradient vector field and performs non-maximum suppression and closed connectivity verification, and outputs a set of lesion voxel boundaries. The mesh reconstruction and consistency verification module performs three-dimensional voxel aggregation, morphological closing operation and surface triangulation reconstruction on the lesion voxel boundary set to obtain a three-dimensional lesion surface mesh, calculates the energy projection value sequence and leakage consistency index and trims and repairs it, and outputs a three-dimensional lesion model.

Images