A method and system for denoising coronary artery calcification score and a storage medium

Abstract

The application provides a coronary artery calcification score denoising method and system and a storage medium, and relates to the technical field of cardiovascular medical image quantitative analysis. The method comprises the following steps: acquiring images of a subject at multiple phases in a resting period; performing spatial positioning without registration on the images at each phase; generating a coronary artery ROI on the image at an anchor phase; labeling high-density voxels in the images at each phase to generate a candidate calcification mask at each phase; performing morphological dilation on the candidate calcification mask at a reference phase; performing an or operation on the candidate calcification mask at the reference phase after the dilation processing on a voxel-by-voxel basis to obtain a confirmation mask; performing an and operation on the candidate calcification mask at the anchor phase without dilation and the confirmation mask on a voxel-by-voxel basis to obtain a final calcification mask; and performing calcification score calculation on the image at the anchor phase based on the final calcification mask. The scheme can reduce the complexity of the process and the deployment threshold.

Description

Technical Field

[0001] This invention relates to the field of quantitative analysis technology of cardiovascular medical images, and in particular to a method, system and storage medium for reducing noise by coronary artery calcium integration. Background Technology

[0002] Coronary Calcification Score (CACS) is a clinically recognized indicator for stratifying the risk of coronary artery disease. The current Agatston score is usually calculated based on 3mm slice thickness images from non-contrast cardiac CT scans. Specifically, it identifies calcifications within the coronary artery pathway with a density not lower than the threshold of 130 HU (Huntersfield units), and accumulates the scores by weighting the lesion area with the peak HU (density). The scores are then used to classify the examinee into five risk categories (0 points / 1–10 points / 11–100 points / 101–400 points / >400 points).

[0003] PCCT (Photon Counting CT) can provide reconstructed images with ultrathin slices of 0.4 mm and others less than 1 mm. Compared with a slice thickness of 3 mm, it has higher spatial resolution and is theoretically more conducive to detecting microcalcifications, thereby improving the ability to display and quantitatively analyze microcalcification candidates.

[0004] However, thinning the slice increases quantum noise. For example, with a slice thickness of 0.4 mm, a large number of random high-density noise voxels appear near the 130 HU threshold, increasing the probability of false positives. Furthermore, the Agatston integral, originally designed for a 3 mm slice thickness, often exhibits systematic bias when directly applied to 0.4 mm ultrathin slice images due to the increased number of slices, the detection of new small lesions, and changes in area statistics.

[0005] Therefore, various solutions have emerged in the industry, such as: A method for averaging or fusing multi-temporal images after registration. This method first performs rigid or non-rigid registration on multiple temporal phases, and then performs averaging, weighted averaging or other fusion to reduce random noise. Deep learning-based multi-temporal denoising methods take multiple temporal phases as network inputs and output denoised images or calcification candidates through neural networks. Such schemes usually require training data, model deployment environment or additional inference resources. Denoising methods at the reconstruction parameter level, such as adjusting layer thickness, reconstruction kernel, virtual single-energy image level, or QIR (quantum iterative reconstruction) intensity, reduce noise in thin-layer images, but do not directly utilize low-correlation-related random noise across time phases.

[0006] While these solutions can be effective to some extent, they still have the following shortcomings: Because the coronary artery diameter is small (about 2-4 mm), the residual error of non-rigid registration under ultra-thin layer conditions can approach or even exceed the lesion scale, resulting in blurred boundaries of microcalcifications, a decrease in HU peak (density), or even falling below the 130 HU threshold. Registration fusion and deep learning inference typically require additional registration optimization, model deployment, or inference environment, and their process complexity is higher than that of direct morphological and logical operation chains. Existing multi-temporal fusion follows the paradigm of (alignment first, averaging later), which has the technical inertia that multi-temporal utilization usually requires image content registration first. There is a lack of a transparent method that can directly utilize the low correlation characteristics of random noise between temporal phases without the need for image content registration. While some deep learning solutions can achieve good noise reduction results, they suffer from problems such as model dependence, weak interpretability, and high deployment environment requirements. Summary of the Invention

[0007] The purpose of this invention is to provide a method, system, and storage medium for reducing noise in coronary artery calcium integration, so as to solve at least one of the above-mentioned technical problems existing in the prior art.

[0008] To facilitate understanding of this invention, before detailing the technical solution, it should be noted that the core of this solution is not averaging after alignment, but rather asymmetric voting that decouples anchor point fidelity from spatial tolerance verification. Its principle includes: Principle 1: PCCT detectors directly count incident X-ray photons, and the noise of each voxel follows a Poisson statistical distribution. In different reconstruction phases of ECG (electrocardiogram) gated acquisition, the random high-density noise at the same spatial location mainly comes from photon count fluctuations and random noise components within the reconstruction window. For different reconstruction phases, the aforementioned random high-density noise can be approximated as independent or low-correlation events.

[0009] Principle 2: Within the cardiac quiescent period determined by ECG signals (e.g., 72%–78% of the RR interval at low heart rates and 30%–40% of the RR interval at high heart rates), adjacent temporal coronary artery displacements can be covered by millimeter-level spatial tolerances. Therefore, it is not necessary to rely on rigid or non-rigid image registration based on image content; cross-temporal verification with limited morphological tolerances is sufficient.

[0010] Principle 3: The anchor point phase is not dilated (to maintain the original spatial accuracy), but only the reference phase is dilated with displacement tolerance (to verify candidate calcifications in the anchor point phase) in an asymmetric structure. This structure avoids morphological expansion of the final scored image itself, and can reduce the impact on the true lesion boundary, HU peak (density), and (calcified) lesion area.

[0011] In a first aspect, to solve the above-mentioned technical problems, the present invention provides a method for reducing noise based on coronary artery calcium integration, comprising the following steps: Step 1: Acquire images of the subject at multiple time points during the quiescent period; the images may be reconstructed images of PCCT ultrathin slices; designate one of the time points as the anchor time point and the remaining time points as reference time points.

[0012] In one feasible implementation, the specific criteria for the quiescent period include at least one of the following criteria: Criterion A: When heart rate < 75 bpm, the mid-to-late diastolic phase is taken as the resting period, such as 72%, 75%, and 78% of the RR interval; when heart rate ≥ 75 bpm, the end-systolic phase is taken as the resting period, such as 30%, 35%, and 40% of the RR interval. Criterion B: The cardiac cycle in which the rate of change of the cross-sectional area (or displacement of the lumen center) of the coronary artery in adjacent time phases is less than a preset threshold is taken as the resting period; the preset threshold is 5% area change (or 0.8 mm displacement). Criterion C: Use the quiescent period specified in the scanning device protocol (or recommended by the reconstruction interface).

[0013] In one feasible implementation, the phase with the most stable cardiac motion, the highest image quality, or the most commonly used phase in clinical practice is selected as the anchor phase.

[0014] Step 2: Based on the patient's physical coordinate system, perform registration-free spatial localization of images at each time phase.

[0015] In one feasible implementation, step 2 includes: Step 21: Read the ImagePositionPatient (patient image position) label from the DICOM (Digital Imaging and Communication in Medicine) image header information of each time phase, and extract the z (patient physical) coordinates of each image layer; at the same time, read spatial geometric information such as ImageOrientationPatient (patient image orientation), PixelSpacing, Rows, Columns, SliceThickness (slice thickness) or inter-slice spacing, and FrameOfReferenceUID (unique identifier of reference coordinate system) to confirm that the images of each time phase are in the same patient physical coordinate system.

[0016] Step 22: For images at each time phase, calculate the common coverage area of ​​the z-coordinate range. The specific calculation formula includes: z_lo = max_i{min(Z_i)}; z_hi = min_i{max(Z_i)}; z_tol = 0.5 × SliceThickness; Where z_lo represents the lower bound of the common coverage area; Z_i represents the set of z coordinates of all slices of the image at the i-th time phase; z_hi represents the upper bound of the common coverage area; z_tol represents the half-slice thickness tolerance; and SliceThickness represents the slice thickness. For images at each time phase, select those that satisfy z_lo The physical mesh is constructed by extracting the corresponding layer sets in order of z_tol≤z≤z_hi+z_tol, and then extracting the layers in order of z coordinate. In this way, we can achieve layer correspondence based on physical space, rather than relying on layer number or image content registration. It should be noted that this step does not perform any resampling, grayscale interpolation, grayscale normalization, grayscale averaging, convolution filtering, or other operations that change the voxel HU values ​​on the original CT grayscale arrays of the images at each time point; this step does not perform rigid registration, affine registration, non-rigid registration, or elastic deformation registration based on the image grayscale content between the images at different time points; the output of this step is the index correspondence between the image layers and voxels in the patient's physical coordinate system for each time point, and this output does not carry any modified grayscale data; Since any resampling or grayscale interpolation of the image array will rewrite the original HU values ​​of the voxels, thereby distorting the calcification segmentation based on the 130 HU threshold and the Agatston density weight based on the HU peak (density), this step is limited to the layer index level, resulting in a physical grid with voxel indexes in the patient's physical coordinate system, so that the original HU values ​​remain unchanged in the cross-temporal spatial correspondence.

[0017] Step 23: If there are differences in the orientation matrix, pixel spacing, matrix size, or field of view of the images at different time points, the physical grid of the image at the anchor point time point is used as the target grid to establish a mapping relationship from the reference time point to the target grid. This will allow the binary candidate calcification mask obtained after the image at the reference time point is filtered by (candidate calcification threshold) to be mapped to the target grid according to the patient's physical coordinates.

[0018] In one feasible implementation, the mapping relationship in step 23 employs nearest neighbor interpolation or an equivalent binary mask resampling method, solely to ensure that subsequent voxel-by-voxel logical operations have the same spatial index basis. This mapping operates on the binary candidate mask rather than the original CT grayscale image; since the mapping object is binary data (with values ​​only 0 or 1), the nearest neighbor mapping does not introduce grayscale interpolation, nor does it modify the original CT image at the anchor point phase. Therefore, the grayscale value used by the final Agatston integral still comes from the image at the anchor point phase (the original reconstructed image).

[0019] Step 24: When the common coverage area is less than a preset range threshold, individual images have missing layers, or there is a (mild) inconsistency in the number of layers, execute any of the following fallback strategies: Backtracking strategy 1: Based on the complete z-coordinate range of the image at the anchor point, images at other time points are matched according to the closest physical plane; Backoff strategy 2: Take the minimum number of common layers with overlapping regions in all time phase images; Backoff strategy 3: For edge (small number) missing layers, remove them only during the final integration, so as not to affect the subsequent voting on the overlapping regions in the middle of the layer; Rollback Strategy 4: When the image of any reference time phase cannot be reliably mapped to the physical grid of the anchor time phase, discard the image of that reference time phase (or output a low confidence warning message).

[0020] Step 3: (In order to exclude non-coronary high-density structures such as vertebral bodies, ribs, aortic walls, and cardiac blood pools) Generate a coronary artery ROI (region of interest) on the image at the anchor point time phase to limit the spatial range for candidate calcification extraction and final integral calculation in subsequent time phase images.

[0021] In one feasible implementation, the mathematical expression for the coronary artery ROI can be: ROI_coronary=Dilate(Mask_myo,R_expand)\(Mask_chamber∪Mask_aorta); Wherein, ROI_coronary represents the coronary artery ROI; Mask_myo represents the myocardial mask; R_expand represents the ROI expansion radius (used to generate the coronary artery candidate region, i.e., the region after expanding the myocardial mask and removing the cardiac blood pools and aorta, consistent across all time phases); Mask_chamber represents the union of the four cardiac blood pool masks; Mask_aorta represents the proximal aortic mask; \ represents the set difference operation; Dilate represents the discrete three-dimensional morphological dilation operation.

[0022] In one feasible implementation, step 3 can be achieved in any of the following ways: Generation Method 1: Using a three-dimensional anatomical segmentation model (e.g., a three-dimensional convolutional network, a Transformer segmentation network, or a rule-model hybrid segmentation algorithm) known in the art, the myocardium, four heart chambers, and the aorta are segmented from the image at the anchor point. The coronary artery ROI is generated according to the aforementioned mathematical expression of the coronary artery ROI. The preferred range for the outer radius of the ROI is 5~10mm. Generation Method 2: First, initial segmentation is performed in the chest transverse image based on threshold ranges for air, soft tissue, and high-density bony structures to remove lung fields, vertebral bodies, ribs, and areas outside the chest wall. Then, within the mediastinal soft tissue connected region, the main cardiac region surrounding the epicardium is determined by combining the image center location, the continuity of the outer edge of the heart, the location of the circular high-density / soft tissue structure near the aorta, and the continuity of the plane. Subsequently, the main cardiac region is expanded outward to form a candidate band-shaped region covering the course of the coronary arteries. Finally, obvious blood pools, aortic walls, valve areas, and bony high-density areas are removed by using HU thresholds, connected region sizes, relative anatomical locations, and plane continuity to obtain the coronary artery ROI. The thresholds and connected region sizes in this method can be set according to pixel spacing, scanning field of view, and local pre-experimentation. When the ROI coverage is insufficient, a low confidence warning message is output or the method is switched to Generation Method 1. Generation Method 3: Simultaneously execute Generation Method 1 and Generation Method 2 to obtain two candidate coronary artery ROIs, namely ROI_A and ROI_B, and obtain the final coronary artery ROI, namely ROI_final, according to any of the following merging strategies: Merging strategy 1: Take the union, i.e., ROI_final=ROI_A∪ROI_B, which is used to supplement the coverage by another method when the coverage of one method is insufficient; Merging strategy 2, taking the intersection, i.e., ROI_final=ROI_A∩ROI_B, is used to perform subsequent voting within the high confidence region agreed upon by both methods; Merging strategy 3: Merge according to layer priority. For example, use ROI_A on layers that can stably provide anatomical segmentation results, fall back to ROI_B on layers that fail to segment or have insufficient coverage, and ensure coverage continuity at cross-layer boundaries through morphological closing operations. It should be noted that the purpose of generation method 3 is to maintain sufficient coverage of the coronary artery course even under conditions where the anatomical segmentation model deployment environment is limited, segmentation at individual levels fails, or image quality distribution is uneven, and to reduce the frequency of low-confidence prompts caused by insufficient coverage of a single method (path). Therefore, the core voting and integration chain of this invention does not depend on a specific deep learning model for implementation.

[0023] In one feasible implementation, the formula for calculating the number of expansion iterations in generation method 1 and generation method 3 is as follows: N_expand=max(1,round(R_expand / PixelSpacing)); Where N_expand represents the number of expansion iterations; PixelSpacing represents the axial pixel spacing; and round represents the rounding operation.

[0024] Step 4: Within the coronary artery ROI, high-density voxels in the images of each time phase are labeled to generate (binarized) candidate calcification masks for each time phase; Specifically, this includes the (binarized) candidate calcification mask Mask_anchor for the anchor phase and the (binarized) candidate calcification mask Mask_ref for the reference phase; The high-density voxels refer to voxels with a density value not lower than a preset density threshold T.

[0025] In one feasible implementation, the preset density threshold T is in the range of 130~200HU.

[0026] Preferably, the preset density threshold T is 130 HU.

[0027] Step 5: Perform morphological dilation on the (binarized) candidate calcification mask of the reference phase.

[0028] In one feasible implementation, step 5 specifically includes: Step 51: Perform three-dimensional morphological dilation only on the (binarized) candidate calcification mask of the reference phase (i.e., do not perform dilation on the candidate calcification mask of the anchor phase) to compensate for small spatial displacements during the static period.

[0029] Step 52: Determine the displacement tolerance D_tol based on the maximum possible physiological displacement D_max (in mm) between the image of the reference phase_i and the image of the anchor point phase; In this way, the physical effect of different relative maximum displacements at different times can be reflected, and the morphological tolerance adapts accordingly. For example, when the anchor point phase is 75% of the RR interval, and the reference phases are 72% and 78%, since the phase difference between these two reference phases and the anchor point phase is 3%, their corresponding D_max may be approximately equal. However, if the reference phase also includes 70% of the interval, then its D_max relative to the anchor point should be determined based on the actual phase difference and the displacement curve during the stationary period, and may be slightly larger than the corresponding values ​​at 72% and 78%.

[0030] The value of D_max shall be derived from at least one of the following sources: The value is derived from: 1. The maximum permissible coronary artery center displacement defined in the quiescent period of the scanning device protocol; Values ​​are derived from: 2. (Patient-specific) RR interval statistics, heart rate fluctuation amplitude, and resting period duration; Value source 3: (quantified in imaging studies) the coronary artery center displacement between images of adjacent cardiac phases; Value source 4: The quiescent period is defined as the change in the cross-sectional area of ​​the coronary artery between adjacent time phases being less than a preset change threshold or the displacement of the coronary artery center being less than 0.8 mm. The maximum allowable displacement of the coronary artery center during this quiescent period is defined as follows: 5. When D_max cannot be directly given by the scanning device protocol (or clinical imaging study), it can be obtained through the inter-phase matching error statistics in the local pre-experiment. Specifically, this includes: identifying a set of cross-phase candidate regions that are clearly the same calcified lesion in the local reference queue; measuring the center displacement of the same calcified lesion between the corresponding images of the anchor phase and the reference phase; and taking the 95% (or 99%) upper quantile of the distribution of the center displacement as the local D_max estimate.

[0031] Thus, based on the above-mentioned value sources, the value range of D_max(phase_i,anchor) is from sub-millimeter to millimeter, with typical values ​​between 0.5 and 2.0 mm; accordingly, the feasible range of D_tol corresponds to 0.5 to 2.0 mm. It should be noted that D_tol is not directly fixed as a single constant based on experience in this field, and can be dynamically adjusted based on any of the following settings: Based on 1, D_max; Based on criterion 2, the subject's heart rate zone; Setting criteria 3: actual spatial resolution and axial pixel spacing of the device; Based on the following criteria: 4. Statistical analysis of the interphase matching error obtained in the preliminary experiment.

[0032] Step 53: Calculate the number of discrete expansion iterations N_tol. The specific formula includes: N_tol=max(1,round(D_tol / PixelSpacing)); PixelSpacing represents the axial pixel spacing; the dilation structure element can be a three-dimensional 6-connected or 26-connected structure element; Step 6: Perform a voxel-by-voxel OR operation (vote counting) on ​​the (all) dilated (processed) candidate calcification masks of the reference time phase (as virtual voters) to obtain the confirmation mask; then perform a voxel-by-voxel AND operation (decision) on the undilated anchor point time phase (as virtual nominators) with the confirmation mask to obtain the final calcification mask; thus, based on the OR / AND voting method, a voxel is included in the final calcification mask only when it is identified as a calcification candidate in the anchor point time phase, and when there is a signal exceeding the threshold (i.e., at least one vote) in the D_tol neighborhood of at least one reference time phase, rather than a majority vote.

[0033] In one feasible implementation, step 6 confirms the specific calculation formula for the mask, including: Mask_confirm= _{i=1..N 1}Dilate(Mask_ref_i,D_tol); Where Mask_confirm represents the confirmation mask; Mask_ref_i represents the (binarized) candidate calcification mask for the i-th reference phase; N represents the total number of phases; This represents finding the union of sets; The specific calculation formula for the final calcification mask is as follows: Mask_final=Mask_anchor∩Mask_confirm; Mask_final represents the final calcification mask.

[0034] In one feasible implementation, step 6 further includes a method for analyzing the inclusion rate of the final calcification mask: Under the condition that random high-density noise at different reference time phases is approximately independent or low-correlated, when N time phases are subjected to asymmetric voting, the probability that any voxel is included in the final calcification mask can be expressed by the upper bound estimation formula: P_FP^{vote}≤p·(1 (1 p_d)^{N 1}); Where p represents the probability that the noise of a single phase exceeds the threshold; p_d represents the probability that the noise of the reference phase accidentally exceeds the threshold again in the neighborhood of D_tol; For example: when p=0.01, p_d=0.02, and N=3, P_FP^{vote}≤0.01× (1 0.98²)≈3.96×10 -4 This demonstrates that cross-temporal confirmation can reduce the probability of isolated, random, high-density candidates being retained; of course, this can also be verified by the number of candidates, integral changes, and risk grading indicators.

[0035] In one feasible implementation, when only two time phases are available (i.e., one anchor time phase and one reference time phase), the specific calculation formula for the final calcification mask degenerates to: Mask_final=Mask_anchor∩Dilate(Mask_ref,D_tol); Although the denoising capability at this time is weaker than that of the three-phase or multi-phase implementation, the low correlation characteristics of cross-phase random noise can still be utilized to reduce isolated random high-density candidates.

[0036] Step 7: Based on the final calcification mask, perform calcification integral calculation on the image at the anchor point.

[0037] In one feasible implementation, step 7 specifically includes: Step 71: Suppose the Lth lesion is composed of a two-dimensional 8-connected component C_L from Mask_final, calculate the lesion area; the specific calculation formula can be expressed as: Area_L=PixelSpacing_x·PixelSpacing_y·#C_L; Where Area_L represents the lesion area; PixelSpacing_x represents the projection of the axial pixel spacing on the x-axis, PixelSpacing_y represents the projection of the axial pixel spacing on the y-axis, and their product is the physical area of ​​a single pixel; #C_L represents the number of pixels contained in C_L.

[0038] Step 72: Let h_L be the peak value of the image pixel HU in the time phase of the anchor point within C_L. The specific calculation formula is: h_L=max_{(x,y)∈C_L}HU_anchor(x,y); The Agatston density weights can be calculated using the following formula: w(h)=1·1{130≤h<200}+2·1{200≤h<300}+3·1{300≤h<400}+4·1{h≥400}; Where w(h) represents the density weight value; 1{·} represents the indicator function.

[0039] Step 73: Filter connected regions with lesion areas smaller than the preset minimum area A_min (e.g., 1.0~1.5 mm²); combine the coronary artery ROI (and conventional relative anatomical location and conventional planar continuity rules), eliminate high-density connected regions located in non-coronary artery anatomical regions (e.g., valve areas, aortic walls, bony structures, etc.), and calculate the (whole heart) total corrected Agatston integral. The specific formula can be expressed as: S_total=Σ_{L∈ _valid}Score_L; Where S_total represents the (whole) total corrected Agatston integral value; _valid represents the set of lesions; Score_L represents the corrected Agatston integral value of the Lth lesion, calculated using the following formula: Score_L=Area_L·w(h_L)·F_thick; Where w(h_L) represents the density weight value of the Lth lesion; F_thick represents the thickness correction factor.

[0040] In one feasible implementation, the setting strategy for F_thick can be any of the following: Strategy 1: Set the theoretical starting point value of F_thick to the ratio of the actual layer thickness to 3 (mm); for example, with a layer thickness of 0.4mm, the theoretical starting point value is 0.133. Strategy 2: (Since ultrathin slice images increase the number of slices and may detect microcalcifications not shown in conventional 3 mm images, directly applying the theoretical starting value is often insufficient to complete the engineering docking with the 3 mm slice thickness Agatston system. Therefore, when a local reference phase queue is available,) F_thick is used as a local calibration parameter and determined by minimizing the error of the clinical 3 mm slice thickness Agatston integral, the consistency of risk grading, or a weighted target of both. For example, a search can be performed within the candidate range to make median(|S_corrected) S_3mm|) or mean(|S_corrected) The smaller value of S_3mm|) is used, that is, the consistency of risk classification is used as a parallel indicator or constraint indicator to minimize the error; where S_corrected represents the corrected whole Agatston integral calculated using F_thick; S_3mm represents the Agatston integral with a 3mm layer thickness; median represents the median; mean represents the arithmetic mean; Strategy 3: Set the implementation range of F_thick to 0.1~0.3.

[0041] Preferably, in the setting strategy 3, the implementation range of F_thick can be set to 0.18~0.22; the specific value can be adjusted according to different devices, reconstruction parameters and local workflows (i.e., combinations of device models, reconstruction parameters, reference queues and data processing flows, etc.).

[0042] It should be noted that the thickness correction method in this scheme differs from the existing correction method that adapts to the standard protocol by adjusting the HU threshold. This correction method retains the 130HU clinical threshold unchanged, and only multiplies the accumulated Score_L by the thickness correction factor F_thick to achieve the magnitude mapping. Therefore, it does not change the extraction rules of candidate voxels near the 130HU threshold, nor does it change the grading boundary of w(h). In contrast, existing methods directly modify the calcification determination boundary at the threshold level (e.g., adjust the threshold from 130HU to 150HU or 170HU), which may change the distribution of candidate voxels near the HU threshold and the determination of the HU peak (density) at the lesion site.

[0043] Step 74: Based on the clinical Agatston risk stratification guidelines, construct five risk level functions, which can be expressed as follows: Risk(s) = {0 if s = 0} 1 if 0 <s≤10 2 if 10 <s≤100 3 if 100 <s≤400 4 if s>400}; Where Risk(s) represents the risk level; s represents S_total.

[0044] In one feasible implementation, the calcification integral denoising method further includes step 8: independently calculating the integral value (i.e., without voting) on ​​the images of each time phase, and calculating the coefficient of variation between time phases.

[0045] In one feasible implementation, the specific formula for calculating the interphase variation coefficient includes: CoV = std({S_i}) / mean({S_i}); Wherein, CoV represents the inter-phase variation coefficient, the lower the coefficient, the better the consistency across phases, and the more reliable the results; std(·) represents the overall standard deviation; S_i represents the single-phase independent Agatston integral of the i-th participating phase; the participating phases include the anchor phase and the reference phase; When mean({S_i}) is less than the preset minimum value If all time-phase integrals are 0, then CoV is recorded as 0; if at least one time-phase integral is not 0, then the CoV value is not calculated, and a low confidence warning message is output. When CoV ≤ 15%, the confidence level is acceptable. When CoV > 30%, the confidence level is deemed unacceptable and requires manual review.

[0046] In one feasible implementation, CoV can also be used as a dynamic quality control signal in any of the following situations: Scenario 1: When CoV exceeds the preset first mutation threshold, D_tol is dynamically adjusted; for example, when CoV is high, D_tol is appropriately increased to improve the consistency between time phases, and steps 5 to 7 are re-executed. Scenario 2: When the CoV exceeds the preset second mutation threshold, the level of the confidence output label is switched; for example, from acceptable to requiring manual review (or not interpretable in a conventional way). Scenario 3: When there is a (significant) difference between candidate anchor phases for CoV, the anchor phase corresponding to the lower value of CoV is used to recalculate the integral result; Scenario 4: In a multi-device protocol or multi-center environment, use the CoV median as one of the objective functions for local parameter calibration (including D_tol, A_min, and F_thick); Scenario 5: When CoV exceeds the third mutation threshold, output an additional label in the structured results indicating whether the CoV state is still supported by the reference phase for each retained lesion.

[0047] Secondly, based on the same inventive concept, this application also provides a coronary artery calcium integration noise reduction system, including a data acquisition module, a data processing module and a result generation module; The data acquisition module is used to acquire images of the subject at multiple time phases during the quiescent period; the images are reconstructed images of PCCT ultrathin slices; one of the time phases is designated as the anchor time phase, and the remaining time phases are reference time phases; The data processing module includes a spatial positioning unit, an ROI generation unit, a mask generation unit, an expansion unit, a mask calculation unit, and an integral calculation unit; The spatial positioning unit is used to perform registration-free spatial positioning of images at each time phase based on the patient's physical coordinate system. The ROI generation unit is used to generate a coronary artery ROI on the image at the anchor point time phase, and to define the spatial range for candidate calcification extraction and final integral calculation of images at subsequent time phases. The mask generation unit is used to mark high-density voxels in images of each time phase within the coronary artery ROI, and generate (binarized) candidate calcification masks for each time phase, including (binarized) candidate calcification mask Mask_anchor for the anchor point phase and (binarized) candidate calcification mask Mask_ref for the reference phase; the high-density voxels refer to voxels with a density value not lower than a preset density threshold T; The expansion unit is used to morphologically expand the (binarized) candidate calcification mask of the reference time phase; The mask operation unit is used to perform a voxel-by-voxel OR operation on the (binarized) candidate calcification masks of the (all) dilated (processed) reference time phases to obtain the confirmed mask; and then perform a voxel-by-voxel AND operation on the (binarized) candidate calcification masks of the undilated anchor point time phases and the confirmed mask to obtain the final calcification mask. The integral calculation unit performs calcification integral calculation on the image at the anchor point time phase based on the final calcification mask; The result generation module is used to send the processing results of the data processing module to external systems.

[0048] In one feasible implementation, the data processing module further includes a coefficient of variation unit, which is used to independently calculate the integral value on the images of each time phase and calculate the coefficient of variation between time phases.

[0049] Thirdly, based on the same inventive concept, this application also provides a computer-readable storage medium, wherein the program stored in the storage medium, when executed by a processor, is able to complete the steps of the aforementioned coronary artery calcium integral noise reduction method.

[0050] By adopting the above technical solution, the present invention has the following beneficial effects: The present invention provides a method, system and storage medium for coronary artery calcium integral noise reduction, which replaces rigid or non-rigid registration based on image grayscale content by using DICOM physical coordinate alignment, binary mask physical mesh mapping and finite displacement tolerance verification, thereby reducing process complexity and deployment threshold. This scheme does not perform any resampling, grayscale interpolation, weighted averaging, or convolution rewriting on the original CT grayscale arrays of the anchor phase and reference phase throughout the entire chain of cross-temporal spatial correspondence, binary mask mapping, reference dilation, logical voting, and final integration; the 130HU threshold and HU peak (density) grading boundary on which the Agatston integral depends are not changed; and the clinical risk grading criteria sensitive to HU peak (density) are preserved in the engineering implementation; This scheme only dilates the image of the reference time phase, avoiding the boundary blurring and HU peak (density) reduction problems caused by morphological expansion or multi-frame averaging of the final integration time phase itself. The parameterization path of this scheme links morphological tolerance to physical displacement and image resolution during the resting period, avoiding empirical constants that are unrelated to the physical resolution of the image.

[0051] This scheme can suppress isolated random high-density candidates without changing the original image density value of the anchor point phase; The core denoising and integration chain of this solution is based on thresholding, masking, dilation, and logical operations, and does not use black-box reasoning results as the voting basis; the coronary artery ROI can be provided by a general anatomical segmentation model, rule-based method, or two-path parallel merging method, so it can be directly deployed on a regular workstation; This scheme uses a thickness correction factor to make the results of ultrathin layers comparable in magnitude to the conventional 3 mm layer thickness Agatston integral risk classification, without changing the 130HU threshold and density weight classification boundary. Based on cross-temporal integral consistency, this scheme uses CoV not only as a result reliability indicator, but also as a dynamic quality control signal within the process to drive displacement tolerance adjustment, anchor point switching, parameter calibration objective function, and reliability level switching. Attached Figure Description

[0052] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0053] Figure 1 A flowchart of a coronary artery calcification integral noise reduction method provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of the common coverage area provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of coronary artery ROI generation provided in an embodiment of the present invention; wherein, diagram a shows generation method 1; and diagram b shows generation method 2. Figure 4 This is a schematic diagram of the final calcification mask generation provided in an embodiment of the present invention. Detailed Implementation

[0054] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0055] In the description of this invention, it should be noted that the terms "center," "upper," "lower," "left," "right," "vertical," "horizontal," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing the invention and for simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on the invention. Furthermore, the terms "first," "second," and "third" are used for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0056] In the description of this invention, it should be noted that, unless otherwise explicitly specified and limited, the terms "installation," "connection," and "linking" should be interpreted broadly. For example, they can refer to a fixed connection, a detachable connection, or an integral connection; they can refer to a mechanical connection or an electrical connection; they can refer to a direct connection or an indirect connection through an intermediate medium; and they can refer to the internal connection of two components. Those skilled in the art can understand the specific meaning of the above terms in this invention based on the specific circumstances.

[0057] The present invention will be further explained below with reference to specific embodiments.

[0058] Example 1: like Figure 1 As shown in the figure, the coronary artery calcium integration noise reduction method provided in this embodiment includes the following steps: Step 1: Acquire images of the subject at multiple time points during the quiescent period; the images are reconstructed images of PCCT ultrathin slices; designate one of the time points as the anchor time point, and the remaining time points as reference time points.

[0059] Furthermore, the specific criteria for the quiescent period include at least one of the following criteria: Criterion A: When heart rate < 75 bpm, the mid-to-late diastolic phase is taken as the resting period, such as 72%, 75%, and 78% of the RR interval; when heart rate ≥ 75 bpm, the end-systolic phase is taken as the resting period, such as 30%, 35%, and 40% of the RR interval. Criterion B: The cardiac cycle in which the rate of change of the cross-sectional area (or displacement of the lumen center) of the coronary artery in adjacent time phases is less than a preset threshold is taken as the resting period; the preset threshold is 5% area change (or 0.8 mm displacement). Criterion C: Use the quiescent period specified in the scanning device protocol (or recommended by the reconstruction interface).

[0060] Furthermore, the phase with the most stable cardiac motion, the highest image quality, or the most commonly used clinical practice was selected as the anchor phase.

[0061] Step 2: Based on the patient's physical coordinate system, perform registration-free spatial localization of images at each time phase.

[0062] Further, step 2 includes: Step 21: Read the ImagePositionPatient (patient image position) label from the DICOM (Digital Imaging and Communication in Medicine) image header information of each time phase, and extract the z (patient physical) coordinates of each image layer; at the same time, read spatial geometric information such as ImageOrientationPatient (patient image orientation), PixelSpacing, Rows, Columns, SliceThickness (slice thickness) or inter-slice spacing, and FrameOfReferenceUID (unique identifier of reference coordinate system) to confirm that the images of each time phase are in the same patient physical coordinate system; When multiple time phases have the same FrameOfReferenceUID, the same orientation matrix, pixel pitch difference less than a preset threshold, and a common interval exists in the Z-axis coverage, the system enters the standard processing flow. If the FrameOfReferenceUID is missing, but the image orientation, pixel pitch, layer coordinates, and scan sequence time information all meet the geometric consistency conditions within the same check, this situation can be marked as a computable state, and the geometric verification identifier can be retained in the processing result. If the orientation matrix is ​​flipped, the pixel pitch difference exceeds the threshold, the layer coordinates are not monotonic, the common coverage is less than the preset number of layers, or any reference time phase cannot be mapped to the anchor point physical grid, then voxel-by-voxel logical voting is not performed; instead, abnormal reference time phases are removed or low confidence warning messages are output. Registration-free spatial localization does not mean ignoring spatial relationships, but rather using the patient's physical coordinates and image geometric labels as the basis for cross-temporal spatial correspondence. By setting a geometric consistency gating before voting, it is possible to prevent images from different fields of view, different matrices, or different directions from being erroneously merged, thereby improving the feasibility of this method in real examination data.

[0063] Step 22, as follows Figure 2 As shown, for images at each time phase, the common coverage area of ​​the z-coordinate range is calculated. The specific calculation formula includes: z_lo = max_i{min(Z_i)}; z_hi = min_i{max(Z_i)}; z_tol = 0.5 × SliceThickness; Where z_lo represents the lower bound of the common coverage area; Z_i represents the set of z coordinates of all slices of the image at the i-th time phase; z_hi represents the upper bound of the common coverage area; z_tol represents the half-slice thickness tolerance; and SliceThickness represents the slice thickness. For images at each time phase, select those that satisfy z_lo The physical mesh is constructed by extracting the corresponding layer sets in order of z_tol≤z≤z_hi+z_tol, and then extracting the layers in order of z coordinate. In this way, we can achieve layer correspondence based on physical space, rather than relying on layer number or image content registration. It should be noted that this step does not perform any resampling, grayscale interpolation, grayscale normalization, grayscale averaging, convolution filtering, or other operations that change the voxel HU values ​​on the original CT grayscale arrays of the images at each time point; this step does not perform rigid registration, affine registration, non-rigid registration, or elastic deformation registration based on the image grayscale content between the images at different time points; the output of this step is the index correspondence between the image layers and voxels in the patient's physical coordinate system for each time point, and this output does not carry any modified grayscale data; Since any resampling or grayscale interpolation of the image array will rewrite the original HU values ​​of the voxels, thereby distorting the calcification segmentation based on the 130 HU threshold and the Agatston density weight based on the HU peak (density), this step is limited to the layer index level, resulting in a physical grid with voxel indexes in the patient's physical coordinate system, so that the original HU values ​​remain unchanged in the cross-temporal spatial correspondence.

[0064] Step 23: If there are differences in the orientation matrix, pixel spacing, matrix size, or field of view of the images at different time points, the physical grid of the image at the anchor point time point is used as the target grid to establish a mapping relationship from the reference time point to the target grid. This will allow the binary candidate calcification mask obtained after the image at the reference time point is filtered by (candidate calcification threshold) to be mapped to the target grid according to the patient's physical coordinates.

[0065] Specifically, when the pixel spacing of the reference phase and the anchor phase is the same or approximately the same, the nearest anchor voxel can be directly mapped according to the physical coordinates; when there is a slight difference in the pixel spacing between the two, the physical coordinate index of the anchor grid can be established first, and then the positive voxel of the reference mask can be projected to the nearest anchor voxel; when the reference voxel is mapped outside the boundary of the anchor grid, the voxel does not participate in the subsequent voting (step 6); when multiple reference voxels are mapped to the same anchor voxel, the anchor voxel still remains binary positive and does not accumulate weight; This approach preserves two boundaries: Boundary 1, spatial correspondence depends on DICOM geometric information and binary mask coordinate transformation, rather than rigid or non-rigid registration of the original CT image based on image content; Boundary 2: The final integration is still performed on the original grayscale image of the anchor point phase. The reference phase only provides neighborhood confirmation information and does not provide the final grayscale value. This boundary helps to preserve the peak density of microcalcifications and the original boundary in the anchor point phase.

[0066] Furthermore, the mapping relationship in step 23 employs nearest neighbor interpolation or an equivalent binary mask resampling method, solely to ensure that subsequent voxel-by-voxel logical operations have the same spatial index basis. This mapping applies to the binary candidate calcification mask rather than the original CT grayscale image; since the mapping object is binary data (taking only 0 or 1 values), the nearest neighbor mapping does not introduce grayscale interpolation, nor does it modify the original CT image at the anchor point phase. Therefore, the grayscale value used by the final Agatston integral still comes from the image at the anchor point phase (the original reconstructed image).

[0067] Step 24: When the common coverage area is less than a preset range threshold, individual images have missing layers, or there is a (mild) inconsistency in the number of layers, execute any of the following fallback strategies: Backtracking strategy 1: Based on the complete z-coordinate range of the image at the anchor point, images at other time points are matched according to the closest physical plane; Backoff strategy 2: Take the minimum number of common layers with overlapping regions in all time phase images; Backoff strategy 3: For layers with a small number of missing edges, remove them only during the final integration, so as not to affect the subsequent voting on the overlapping regions in the middle of the layer; Rollback Strategy 4: When the image of any reference time phase cannot be reliably mapped to the physical grid of the anchor time phase, discard the image of that reference time phase (or output a low confidence warning message).

[0068] Step 3: (In order to exclude non-coronary high-density structures such as vertebral bodies, ribs, aortic walls, and cardiac blood pools) Generate a coronary artery ROI (region of interest) on the image at the anchor point time phase to limit the spatial range for candidate calcification extraction and final integral calculation in subsequent time phase images.

[0069] Furthermore, the mathematical expression for the coronary artery ROI can be: ROI_coronary=Dilate(Mask_myo,R_expand)\(Mask_chamber∪Mask_aorta); Wherein, ROI_coronary represents the coronary artery ROI; Mask_myo represents the myocardial mask; R_expand represents the ROI expansion radius; Mask_chamber represents the union of the four cardiac chamber blood pool masks; Mask_aorta represents the proximal aortic mask; \ represents the set difference operation; and Dilate represents the discrete three-dimensional morphological dilation operation.

[0070] Furthermore, such as Figure 3 As shown, step 3 can be implemented using any of the following generation methods: Generation method 1 (anatomical segmentation model), such as Figure 3 As shown in Figure a, the myocardium, four cardiac chambers, and aorta are segmented from the image at the anchor point using a three-dimensional anatomical segmentation model (such as a three-dimensional convolutional network, a Transformer segmentation network, or a rule-model hybrid segmentation algorithm) known in the art. The coronary artery ROI is generated according to the aforementioned mathematical expression of the coronary artery ROI. The preferred range of the ROI outer radius is 5~10mm. Generation method 2 (rule-based), such as Figure 3As shown in Figure b, initial segmentation is first performed in the chest transverse image based on threshold ranges for air, soft tissue, and high-density bony structures to remove the lung field, vertebral body, ribs, and extrathoracic wall regions. Then, within the mediastinal soft tissue connected region, the main cardiac region surrounding the epicardium is determined by combining the image center location, the continuity of the outer edge of the heart, the location of the circular high-density / soft tissue structure near the aorta, and the layer continuity. Subsequently, the main cardiac region is expanded outward to form a candidate band-shaped region covering the course of the coronary arteries. Finally, obvious blood pools, aortic walls, valve areas, and bony high-density areas are removed by using HU thresholds, connected region sizes, relative anatomical locations, and layer continuity to obtain the coronary artery ROI. The thresholds and connected region sizes in this method can be set according to pixel spacing, scanning field of view, and local pre-experimentation. When the ROI coverage is insufficient, a low confidence warning message is output or the generation method 1 is switched. Generation Method 3: Simultaneously execute Generation Method 1 and Generation Method 2 to obtain two candidate coronary artery ROIs, namely ROI_A and ROI_B, and obtain the final coronary artery ROI, namely ROI_final, according to any of the following merging strategies: Merging strategy 1: Take the union, i.e., ROI_final=ROI_A∪ROI_B, which is used to supplement the coverage by another method when the coverage of one method is insufficient; Merging strategy 2, taking the intersection, i.e., ROI_final=ROI_A∩ROI_B, is used to perform subsequent voting within the high confidence region agreed upon by both methods; Merging strategy 3: Merge according to layer priority. For example, use ROI_A on layers that can stably provide anatomical segmentation results, fall back to ROI_B on layers that fail to segment or have insufficient coverage, and ensure coverage continuity at cross-layer boundaries through morphological closing operations. The criteria for prioritizing levels include: Evaluation criteria 1: Whether the anatomical segmentation model outputs a complete myocardial mask at this level; Assessment Criterion 2: Does the rule-based method output a continuous outer edge of the heart at this level? Assessment criterion 3: Whether this level falls within typical boundary areas such as the apex, base, or root of the heart; Evaluation Criterion 4: Whether the continuity of ROI between this layer and adjacent layers satisfies the spatial consistency constraint.

[0071] It should be noted that the purpose of generation method 3 is to maintain sufficient coverage of the coronary artery course even under conditions where the anatomical segmentation model deployment environment is limited, segmentation at individual levels fails, or image quality distribution is uneven, and to reduce the frequency of low-confidence prompts caused by insufficient coverage of a single method (path). Therefore, the core voting and integration chain of this invention does not depend on a specific deep learning model for implementation.

[0072] Furthermore, in generation methods 1 and 3, the formula for calculating the number of expansion iterations is as follows: N_expand=max(1,round(R_expand / PixelSpacing)); Where N_expand represents the number of expansion iterations; round represents the rounding operation.

[0073] Step 4: Within the coronary artery ROI, high-density voxels in the images of each time phase are marked to generate binary candidate calcification masks for each time phase. Specifically, this includes the binary candidate calcification mask Mask_anchor for the anchor phase and the binary candidate calcification mask Mask_ref for the reference phase; The high-density voxels refer to voxels with a density value not lower than a preset density threshold T.

[0074] Furthermore, the preset density threshold T ranges from 130 to 200 HU.

[0075] Preferably, the preset density threshold T is 130 HU.

[0076] Step 5: Perform morphological dilation on the binarized candidate calcification mask of the reference phase.

[0077] Furthermore, step 5 specifically includes: Step 51: Perform three-dimensional morphological dilation only on the binarized candidate calcification mask of the reference time phase (i.e., do not perform dilation on the binarized candidate calcification mask of the anchor point time phase) to compensate for small spatial displacements during the static period.

[0078] Step 52: Determine the displacement tolerance D_tol based on the maximum possible physiological displacement D_max (in mm) between the image of the reference phase_i and the image of the anchor point phase; In this way, the physical effect of different relative maximum displacements at different times can be reflected, and the morphological tolerance adapts accordingly. For example, when the anchor point phase is 75% of the RR interval, and the reference phases are 72% and 78%, since the phase difference between these two reference phases and the anchor point phase is 3%, their corresponding D_max may be approximately equal. However, if the reference phase also includes 70% of the interval, then its D_max relative to the anchor point should be determined based on the actual phase difference and the displacement curve during the stationary period, and may be slightly larger than the corresponding values ​​at 72% and 78%.

[0079] The value of D_max shall be derived from at least one of the following sources: The value is derived from: 1. The maximum permissible coronary artery center displacement defined in the quiescent period of the scanning device protocol; Values ​​are derived from: 2. (Patient-specific) RR interval statistics, heart rate fluctuation amplitude, and resting period duration; Value source 3: (quantified in imaging studies) the coronary artery center displacement between images of adjacent cardiac phases; Value source 4: The quiescent period is defined as the change in the cross-sectional area of ​​the coronary artery between adjacent time phases being less than a preset change threshold or the displacement of the coronary artery center being less than 0.8 mm. The maximum allowable displacement of the coronary artery center during this quiescent period is defined as follows: 5. When D_max cannot be directly given by the scanning device protocol (or clinical imaging study), it can be obtained through the inter-phase matching error statistics in the local pre-experiment. Specifically, this includes: identifying a set of cross-phase candidate regions that are clearly the same calcified lesion in the local reference queue; measuring the center displacement of the same calcified lesion between the corresponding images of the anchor phase and the reference phase; and taking the 95% (or 99%) upper quantile of the distribution of the center displacement as the local D_max estimate.

[0080] Thus, based on the above-mentioned value sources, the value range of D_max(phase_i,anchor) is from sub-millimeter to millimeter, with typical values ​​between 0.5 and 2.0 mm; accordingly, the feasible range of D_tol is set to 0.5 to 2.0 mm. It should be noted that D_tol is not directly fixed as a single constant based on experience in this field, and can be dynamically adjusted based on any of the following settings: Based on 1, D_max; Based on criterion 2, the subject's heart rate zone; Setting criteria 3: actual spatial resolution and axial pixel spacing of the device; Based on the following criteria: 4. Statistical analysis of inter-phase matching errors obtained in the preliminary experiment; Step 53: Calculate the number of discrete expansion iterations N_tol. The specific calculation formula includes: N_tol=max(1,round(D_tol / PixelSpacing)); PixelSpacing represents the axial pixel spacing; the dilation structure element can be a three-dimensional 6-connected or 26-connected structure element; Step 6: Perform a voxel-by-voxel OR operation (vote counting) on ​​all binarized candidate calcification masks of the reference time phases after dilation (as virtual voters) to obtain the confirmed mask; then perform a voxel-by-voxel AND operation (decision) on the binarized candidate calcification masks of the anchor time phases (as virtual nominators) with the confirmed mask to obtain the final calcification mask; thus, based on the OR / AND voting method, a voxel is included in the final calcification mask only if it is identified as a calcification candidate in the anchor time phase and there is a signal exceeding the threshold (i.e., at least one vote) in the D_tol neighborhood of at least one reference time phase, rather than a majority vote.

[0081] Furthermore, such as Figure 4 As shown, in step 6, the specific calculation formula for the mask is confirmed to include: Mask_confirm= _{i=1..N 1}Dilate(Mask_ref_i,D_tol); Where Mask_confirm represents the confirmation mask; Mask_ref_i represents the binary candidate calcification mask for the i-th reference phase; N represents the total number of phases; The specific calculation formula for the final calcification mask is as follows: Mask_final=Mask_anchor∩Mask_confirm; Mask_final represents the final calcification mask.

[0082] Furthermore, step 6 also includes a method for analyzing the inclusion rate of the final calcification mask: Under the condition that random high-density noise at different reference time phases is approximately independent or low-correlated, when N time phases are subjected to asymmetric voting, the probability that any voxel is included in the final calcification mask can be expressed by the upper bound estimation formula: P_FP^{vote}≤p·(1 (1 p_d)^{N 1}); Where p represents the probability that the noise of a single phase exceeds the threshold; p_d represents the probability that the noise of the reference phase accidentally exceeds the threshold again in the neighborhood of D_tol; For example: when p=0.01, p_d=0.02, and N=3, P_FP^{vote}≤0.01× (1 0.98²)≈3.96×10 -4 This demonstrates that cross-temporal confirmation can reduce the probability of isolated, random, high-density candidates being retained; of course, this can also be verified by the number of candidates, integral changes, and risk grading indicators.

[0083] Furthermore, when only two time phases are available (i.e., one anchor time phase and one reference time phase), the specific calculation formula for the final calcification mask degenerates into: Mask_final=Mask_anchor∩Dilate(Mask_ref,D_tol); Although the denoising capability at this time is weaker than that of the three-phase or multi-phase implementation, the low correlation characteristics of cross-phase random noise can still be utilized to reduce isolated random high-density candidates.

[0084] Step 7: Based on the final calcification mask, perform calcification integral calculation on the image at the anchor point.

[0085] Furthermore, step 7 specifically includes: Step 71: Suppose the Lth lesion is composed of a two-dimensional 8-connected component C_L from Mask_final, calculate the lesion area; the specific calculation formula can be expressed as: Area_L=PixelSpacing_x·PixelSpacing_y·#C_L; Where Area_L represents the lesion area; PixelSpacing_x represents the projection of the axial pixel spacing on the x-axis, PixelSpacing_y represents the projection of the axial pixel spacing on the y-axis, and their product is the physical area of ​​a single pixel; #C_L represents the number of pixels contained in C_L.

[0086] Step 72: Let h_L be the peak value of the image pixel HU in the time phase of the anchor point within C_L. The specific calculation formula is: h_L=max_{(x,y)∈C_L}HU_anchor(x,y); The Agatston density weights can be calculated using the following formula: w(h)=1·1{130≤h<200}+2·1{200≤h<300}+3·1{300≤h<400}+4·1{h≥400}; Where w(h) represents the density weight value; 1{·} represents the indicator function.

[0087] Step 73: Filter connected regions with lesion areas smaller than the preset minimum area A_min (e.g., 1.0~1.5 mm²); combine the coronary artery ROI (and conventional relative anatomical location and conventional planar continuity rules), remove high-density connected regions located in non-coronary anatomical regions (e.g., valve areas, aortic walls, bony structures, etc.), and calculate the (whole heart) total corrected Agatston integral. The specific calculation formula can be expressed as: S_total=Σ_{L∈ _valid}Score_L; Where S_total represents the (whole) total corrected Agatston integral value; _valid represents the set of lesions; Score_L represents the corrected Agatston integral value of the Lth lesion, calculated using the following formula: Score_L=Area_L·w(h_L)·F_thick; Where w(h_L) represents the density weight value of the Lth lesion; F_thick represents the thickness correction factor.

[0088] Furthermore, the setting strategy for F_thick can be any of the following: Strategy 1: Set the theoretical starting point value of F_thick to the ratio of the actual layer thickness to 3 (mm); for example, with a layer thickness of 0.4mm, the theoretical starting point value is 0.133. Strategy 2: (Since ultrathin layer images increase the number of layers and may detect microcalcifications not shown in conventional 3 mm images, directly applying the theoretical starting value is often insufficient to complete the engineering docking with the conventional 3 mm thick Agatston system. Therefore, when a local reference phase queue is available,) F_thick is used as a local calibration parameter and determined by minimizing the error of the integral with the conventional 3 mm thick Agatston system, the consistency of risk classification, or a weighted objective of the two. For example, a search can be performed within the candidate range to make median(|S_corrected) S_3mm|) or mean(|S_corrected) The smaller value of S_3mm|) is used, that is, the consistency of risk classification is used as a parallel indicator or constraint indicator to minimize the error; where S_corrected represents the corrected whole Agatston integral calculated using F_thick; S_3mm represents the Agatston integral with a 3mm layer thickness; median represents the median; mean represents the arithmetic mean; Strategy 3: Set the implementation range of F_thick to 0.1~0.3.

[0089] Preferably, in the setting strategy 3, the implementation range of F_thick can be set to 0.18~0.22; the specific value can be adjusted according to different devices, reconstruction parameters and local workflows.

[0090] It should be noted that the thickness correction method in this scheme differs from the existing correction method that adapts to the standard protocol by adjusting the HU threshold. This correction method retains the 130HU clinical threshold unchanged, and only multiplies the accumulated Score_L by the thickness correction factor F_thick to achieve the magnitude mapping. Therefore, it does not change the extraction rules of candidate voxels near the 130HU threshold, nor does it change the grading boundary of w(h). In contrast, existing methods directly modify the calcification determination boundary at the threshold level (e.g., adjust the threshold from 130HU to 150HU or 170HU), which may change the distribution of candidate voxels near the HU threshold and the determination of the HU peak (density) at the lesion site.

[0091] Step 74: Based on the clinical Agatston risk stratification guidelines, construct five risk level functions, which can be expressed as follows: Risk(s) = {0 if s = 0} 1 if 0 <s≤10 2 if 10 <s≤100 3 if 100 <s≤400 4 if s>400}; Here, Risk(s) represents the risk level, which is a numerical range label and is an intermediate quantitative result for doctors' reference, but does not constitute a disease diagnosis conclusion; s represents S_total.

[0092] Furthermore, this calcification integral denoising method also includes step 8: independently calculating the integral value (i.e., without voting) on ​​the images of each time phase, and calculating the coefficient of variation between time phases.

[0093] Furthermore, the specific formula for calculating the interphase variation coefficient includes: CoV = std({S_i}) / mean({S_i}); Wherein, CoV represents the inter-phase variation coefficient, the lower the coefficient, the better the consistency across phases, and the more reliable the results; std(·) represents the overall standard deviation; S_i represents the single-phase independent Agatston integral of the i-th participating phase; the participating phases include the anchor phase and the reference phase; When mean({S_i}) is less than the preset minimum value If all time-phase integrals are 0, then CoV is recorded as 0; if at least one time-phase integral is not 0, then the CoV value is not calculated, and a low confidence warning message is output. When CoV ≤ 15%, the confidence level is acceptable. When CoV > 30%, the confidence level is deemed unacceptable and requires manual review.

[0094] Furthermore, CoV can also serve as a dynamic quality control signal, applicable to any of the following situations: Scenario 1: When CoV exceeds the preset first mutation threshold, D_tol is dynamically adjusted; for example, when CoV is high, D_tol is appropriately increased to improve the consistency between time phases, and steps 5 to 7 are re-executed. If CoV decreases to an acceptable range after recalculation, output the standard result and record the D_tol used in the quality control label; if CoV still does not decrease, do not increase D_tol (to avoid decoupling from the upper limit of cardiac physiological displacement), and proceed to scenario 2 or scenario 3.

[0095] Scenario 2: When the CoV exceeds the preset second mutation threshold, the level of the confidence output label is switched; for example, from acceptable to requiring manual review (or not interpretable in a conventional way). Switching levels does not change the rules for calculating points, but it will trigger graphical interface prompts and corresponding fields in the structured results file. Scenario 3: When there is a (significant) difference between candidate anchor phases for CoV, the anchor phase corresponding to the lower value of CoV is used to recalculate the integral result; Anchor point switching does not average the results of multiple anchor points; it only uses the independent results corresponding to the lower anchor point of CoV. Scenario 4: In a multi-device protocol or multi-center environment, use the CoV median as one of the objective functions for local parameter calibration (including D_tol, A_min, and F_thick); It can be combined with integral error and risk stratification consistency to form a multi-objective weighted system. For example, the calculation formula for calibration objective J can be defined as follows: J(D_tol,A_min,F_thick)=α·median(|S_corrected S_3mm|)+β·(1 concordance(Risk_corrected, Risk_3mm))+γ·median(CoV); Where α, β, and γ are local weights; concordance represents the consistency measure between the two risk grading results (ranging from 0 to 1, with 1 indicating complete consistency); Risk_corrected represents the risk grade obtained by S_corrected according to the risk grade function; and Risk_3mm represents the risk grade obtained by the Agatston integral of (clinical) 3mm slice thickness according to the risk grade function. This calculation formula ensures that parameter selection not only serves the integration error and hierarchical consistency, but also takes into account cross-temporal phase consistency. Scenario 5: When CoV exceeds the third mutation threshold, an additional label is output in the structured results indicating whether each preserved lesion is still supported by the reference phase in this CoV state. This label helps users determine whether the anomalous CoV is driven by a few large-area candidates or local high-density noise clusters.

[0096] Example 2: Based on the method in Example 1, a PCCT non-contrast cardiac scan experiment of a subject with a low heart rate will be briefly described.

[0097] Reconstruction parameters: Layer thickness 0.4mm; QIR: Level 2; Virtual single-energy image (VMI): 70 keV; The mid-to-late diastolic phase was taken as the quiescent phase; Phase selection: 72%, 75%, and 78%, a total of three phases (i.e., three series); 75% of the time phase was selected as the anchor time phase, and 72% and 78% of the time phase were selected as reference time phases; Read the ImagePositionPatient Z coordinates of three DICOM sets, calculate the common coverage area according to the common coverage area calculation formula, and extract the common layer set (assuming 800 layers). On the image at 75% of the time phase, the myocardium and cardiac chamber masks are obtained according to generation method 1, with R_expand=8mm and PixelSpacing=0.39mm; according to the formula for calculating the number of dilation iterations, N_expand=round(8 / 0.39)=21 dilation iterations are calculated. Within the coronary artery ROI, candidate calcification masks Mask_72, Mask_75, and Mask_78 with HU≥130 were extracted for three time phases. Inflation was performed on Mask_72 and Mask_78 with D_tol=1.0mm. According to the formula for calculating the number of discrete inflation iterations, N_tol=3 times was obtained. Mask_75 does not expand; where D_tol=1.0mm is derived from the maximum physiological displacement of the reference phase (72%, 78%) relative to the anchor phase (75%) during the rest period; By confirming the calculation formula for the mask and the specific calculation formula for the final calcified mask, we can calculate: Mask_confirm = Dilated_72 ∪ Dilated_78; Mask_final = Mask_75 ∩ Mask_confirm; Label the 2D 8-connected components of each layer in Mask_final. Two valid lesions were detected in a certain layer: Lesion 1: According to the lesion area calculation formula, the lesion area is 2.1 mm², and the HU peak value is 268 HU. According to the Agatston density weight calculation formula, w=2. According to the modified Agatston integral calculation formula, the contribution is Score_L=2.1×2×0.20=0.84. For lesion 2, similarly, we can obtain: the lesion area is 5.8 mm², the peak HU value is 412 HU, w = 4, and the contribution is 5.8 × 4 × 0.2 = 4.64. Then, according to the formula for calculating the total modified Agatston integral, S_total is obtained, and then according to the five risk level functions, Risk(s) is obtained; The Agatston integrals are calculated separately at the three time phases: S_72 = 45.2, S_75 = 48.1, S_78 = 43.7; According to the formula for calculating the coefficient of variation between time phases, CoV = 1.83 / 45.67 ≈ 4.0%, indicating good consistency between time phases.

[0098] Throughout the entire process of this embodiment, the CT grayscale array of the original anchor phase image was not resampled or grayscale interpolated; the HU peak and lesion area used for the final integration were both derived from the original reconstructed image of the anchor phase (75%).

[0099] Example 3: For subjects with higher heart rates, the resting phase is switched to the end of systole, for example, 30%, 35%, or 40% of the RR interval; other steps are the same as in Example 2, but D_tol can be appropriately relaxed to 1.2~1.5 mm to cover a larger expected displacement. This relaxation is based on the fact that under high heart rate conditions, the maximum physiological displacement of the reference phase relative to the anchor phase in the end of systole may be higher than that in the mid-to-late diastole (at low heart rates). This indicates that this scheme is not limited to the diastolic phase, and the value of D_tol can be dynamically adjusted according to the heart rate interval and the definition of the resting phase.

[0100] Example 4: Based on Example 1, local layer thickness calibration and single-time comparison verification were performed.

[0101] Cardiac scan data were collected from 52 enrolled patients using PCCT with a slice thickness of 0.4 mm. The clinical reference was the Agatston score with a slice thickness of 3 mm on the same machine. The subjects' ages ranged from 31 to 77 years (mean 55.6 years), with 33 males and 19 females. Their heart rates ranged from 47 to 107 bpm (mean 71.4 bpm).

[0102] The parameter range is shown in Table 1: Table 1

[0103] The results of the three-phase asymmetric voting method were compared with the clinical Agatston integral for a 3mm slice thickness. The results are detailed in Table 2. Table 2

[0104] The three-phase asymmetric voting results of this method were compared with the conventional single-phase 0.4mm layer thickness Agatston integral. After voting (step 6), the number of false positives of the candidate high-density voxels of this method decreased significantly, and the number of true positives was basically consistent with that of the conventional single-phase 0.4mm layer thickness Agatston integral.

[0105] After local calibration of the slice thickness (thickness correction factor F_thick), this method showed a complete agreement rate of 40.4% in comparing the corrected Agatston score for 0.4 mm slice thickness with the clinical Agatston score for 3 mm slice thickness in 5 risk categories, and the proportion of cases with a difference of no more than 1 level was 50.0%; the median coefficient of variation (CoV) between time phases was approximately 7.2%.

[0106] This embodiment demonstrates that the method can verify candidate calcifications of anchor point phases by leveraging the low correlation characteristics of multi-temporal random noise without relying on image content registration. This reduces the burden of isolated random high-density candidates in single-temporal ultrathin layer integration while preserving the spatial accuracy of anchor point phases.

[0107] Example 5: Based on Embodiment 1, step 3, generating method 2 (rule-based method), may specifically include: First, in the anchor point phase image, air, soft tissue and high-density bony structures are distinguished according to the HU value range, and the air area outside the chest wall and obvious bony structure area are removed. Subsequently, candidate regions of the heart body are determined in the mediastinal soft tissue connectivity region, which can be jointly defined by the image center location, the inner contour of the thoracic cavity, the continuity of the outer edge of the heart, and the anatomical location of the proximal aorta. After obtaining the candidate region of the heart body, the region is spatially expanded to cover the epicardial fat and the coronary artery course; the expansion radius can be converted into the number of expansion iterations based on the pixel spacing and voxel size; Subsequently, the obvious cardiac blood pools, aortic lumen, high-density areas adjacent to the aortic wall, valve areas, and bony remnants are subtracted from the expanded area. The subtraction method can be completed by combining location regularity, layer continuity, connected domain size, relative distance from the aortic center, relative position from the boundary of the main heart body, and cross-layer continuity features of high-density structures. This method does not require precise coronary artery centerlines or segmentation of every coronary artery lumen. The goal is to generate a candidate space that covers the possible course of the coronary arteries while excluding most non-coronary high-density structures. Subsequent multi-temporal asymmetric voting and Agatston integration are performed within this candidate space. When the ROI coverage obtained by this method is too wide, subsequent filtering based on connected region area, anatomical location, and slice continuity can continue to eliminate non-coronary structures. When the ROI coverage obtained by this method is insufficient, an insufficient ROI coverage warning message can be output, or the method can switch to the anatomical segmentation model (i.e., generation method 1).

[0108] Example 6: This embodiment provides a coronary artery calcium integration noise reduction system, including a data acquisition module, a data processing module, and a result generation module; The data acquisition module is used to acquire images of the subject at rest during the cardiac motion resting period; the images are reconstructed images of PCCT ultra-thin slices; one of the phases is designated as the anchor phase, and the remaining phases are reference phases; The data processing module includes a spatial positioning unit, an ROI generation unit, a mask generation unit, an expansion unit, a mask calculation unit, and an integral calculation unit; The spatial positioning unit is used to perform registration-free spatial positioning of images at each time phase based on the patient's physical coordinate system. The ROI generation unit is used to generate a coronary artery ROI on the image at the anchor point time phase, and to define the spatial range for candidate calcification extraction and final integral calculation of images at subsequent time phases. The mask generation unit is used to mark high-density voxels in images of each time phase within the coronary artery ROI, and generate binary candidate calcification masks for each time phase, including a binary candidate calcification mask Mask_anchor for the anchor point phase and a binary candidate calcification mask Mask_ref for the reference phase; the high-density voxels refer to voxels with a density value not lower than a preset density threshold T; The expansion unit is used to perform morphological expansion on the binarized candidate calcification mask of the reference time phase. The mask operation unit is used to perform a voxel-by-voxel OR operation on all the binarized candidate calcification masks of the reference time phase after dilation to obtain the confirmed mask; and then perform a voxel-by-voxel AND operation on the binarized candidate calcification masks of the undilated anchor point time phase and the confirmed mask to obtain the final calcification mask. The integral calculation unit performs calcification integral calculation on the image at the anchor point time phase based on the final calcification mask; The result generation module is used to send the processing results of the data processing module to external systems.

[0109] Furthermore, the data processing module also includes a coefficient of variation unit, which is used to independently calculate the integral value on the images of each time phase and calculate the coefficient of variation between time phases.

[0110] Example 7: Based on the noise reduction system mentioned in Example 6, the processing result includes the total corrected Agatston integral value and quality control label.

[0111] The quality control labels include: whether the geometric verification is passed, the number of reference phases participating in the voting, the coverage of the common Z-axis, whether the coronary artery ROI is adequately covered, the independent integral value of each phase, the CoV value, whether there are extreme single-phase outliers, heart rate fluctuation prompts, and credibility prompts, etc.

[0112] The quality control labels can be displayed in a graphical interface or written to a structured results file.

[0113] The credibility information includes: When the CoV level is below a preset threshold, the geometric verification passes, and the coronary artery ROI coverage is sufficient, the system outputs a standard confidence level label. When the CoV threshold is exceeded or a certain reference time phase differs too much from other time phases, the system outputs a check mark. When public coverage is insufficient, all reference phases are unavailable, or coronary artery ROI generation fails, the system outputs uninterpretable markers and retains single-phase integrals for user reference.

[0114] The quality control labels do not replace the doctor's diagnosis, but rather provide indications of the stability of the image processing results and the computational conditions.

[0115] The system's graphical interface simultaneously displays the binarized candidate calcification masks for the anchor phase, the reference phase, the confirmation mask, the final calcification mask, and the removed binarized candidate calcification masks. Users can view the ultimately preserved lesions and their spatial neighborhoods supported by the reference phase layer by layer. This interface helps identify which candidate lesions contribute to the final integral (total modified Agatston integral) and also helps detect problems such as insufficient coronary artery ROI coverage, non-coronary high-density residues, or geometric mapping abnormalities.

[0116] Example 8: The noise reduction system mentioned in Example 6 can be deployed as a CT workstation plugin, stand-alone image post-processing software, in-hospital server program, or post-processing module of an image archiving system.

[0117] When the noise reduction system runs as a server program, it can process multiple case samples in batches and save the structured results in CSV, JSON, DICOM SR or other formats readable by the hospital system.

[0118] When the noise reduction system is run as a workstation plugin, it can automatically identify the available time phases after the user selects a multi-time phase sequence and display the masks for each step on the results interface.

[0119] Regardless of the deployment method used, the core processing chain of this invention remains a transparent process consisting of thresholding, binary masking, morphological dilation, and logical operations.

[0120] Example 9: This embodiment provides a computer-readable storage medium. When the program stored in the storage medium is executed by a processor, it can complete the steps of the aforementioned coronary artery calcium integral noise reduction method.

[0121] Example 10: Based on the calcification integral noise reduction method in Example 1, reconstructed images of ultrathin layers less than 1 mm thick, other than 0.4 mm layer thickness, are processed, such as images with layer thicknesses of 0.2 mm, 0.25 mm, 0.5 mm, 0.6 mm, 0.75 mm, or 0.8 mm.

[0122] Based on the difference between SliceThickness and ImagePositionPatient in the DICOM image header information, the actual slice thickness or reconstruction interval is determined, and the common Z-axis coverage, z_tol, F_thick candidate range and lesion area filtering parameters (i.e., A_min) are adjusted accordingly.

[0123] Example 11: Based on Example 7, for boundary case samples where the total modified Agatston score is close to the cutoff values ​​of 0, 10, 100, or 400, the system can simultaneously output the original score and risk level, and indicate whether the result is close to the cutoff value. For example, when the total modified Agatston score is in the vicinity of 9-11, 95-105, or 380-420, the system can output a boundary prompt so that users can view candidate lesions and quality control labels. This prompt only reflects that small numerical changes near the stratification boundary may affect the category, and does not change the score calculation rules.

[0124] Example 12: Based on Example 1, specific examples of discrete expansion iteration count N_tol include: Example 1: The anchor point phase is 75% phase, PixelSpacing = 0.39 mm, and the reference phases are 72% and 78% phases. The phase difference between the two reference phases and the anchor point is 3%. According to the coronary artery center displacement curve during the quiescent period, D_max(72%, 75%) ≈ D_max(78%, 75%) ≈ 1.0 mm. The calculation yields: N_tol = round(1.0 / 0.39) = 3 discrete expansions.

[0125] Example 2: The anchor point phase is 75% phase, PixelSpacing = 0.39 mm, and the reference phases are 70% and 80% phases. The phase difference between the two reference phases and the anchor point is 5%. According to the coronary artery center displacement curve during the quiescent period, D_max(70%,75%)≈D_max(80%,75%)≈ 1.5 mm. The calculation yields: N_tol = round(1.5 / 0.39) = 4 discrete expansions.

[0126] Calculation example 3: The anchor point phase is 75% of the time phase, with PixelSpacing = 0.39 mm. Reference phases include 72%, 78%, and 80% of the time phase. In this example, the phase difference between the reference phase 80% and the anchor point is larger than that of 72% and 78%, so D_max(80%, 75%) > D_max(72%, 75%) = D_max(78%, 75%). Calculations are performed independently for each time phase. N_tol = 3 is applied to the 72% and 78% phases, while N_tol = 4 is applied to the 80% phase. This example demonstrates the physical meaning that the maximum displacement can differ between different time phases, and the morphological tolerance should adapt accordingly.

[0127] Calculation example 4: Scenario implemented during the end-systolic phase at high heart rate. The anchor point phase is the 35% phase, PixelSpacing = 0.40 mm, and the reference phases are the 30% and 40% phases. Because under high heart rate conditions, the maximum physiological displacement of the reference phase relative to the anchor point during the end-systolic resting phase may be higher than that during the mid-to-late diastolic phase at low heart rate, using D_max(30%, 35%) ≈ D_max(40%, 35%) ≈ 1.3 mm, the calculation yields: N_tol = round(1.3 / 0.40) = 3 discrete expansions.

[0128] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features therein. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for noise reduction based on coronary artery calcium integration, characterized in that, include: Step 1: Acquire images of the subject at multiple time points during the resting period; One of the time phases is designated as the anchor time phase, and the remaining time phases are designated as reference time phases; Step 2: Based on the patient's physical coordinate system, perform registration-free spatial localization of images at each time phase; Step 3: Generate a coronary artery ROI on the image at the anchor point time phase; Step 4: Within the coronary artery ROI, high-density voxels in the images of each time phase are marked to generate candidate calcification masks for each time phase; Step 5: Perform morphological dilation on the candidate calcification mask of the reference phase; Step 6: Perform a voxel-by-voxel OR operation on the candidate calcification masks of the expanded reference phase to obtain the confirmed mask; Then, the candidate calcification mask of the unexpanded anchor point phase is ANDed with the confirmed mask on a voxel-by-voxel basis to obtain the final calcification mask; Step 7: Based on the final calcification mask, perform calcification integral calculation on the image at the anchor point.

2. The calcification integral noise reduction method according to claim 1, characterized in that, The specific criteria for the quiescent period include at least one of the following criteria: Criterion A: When the heart rate is <75 bpm, the mid-to-late diastolic phase is taken as the resting phase; when the heart rate is ≥75 bpm, the end-systolic phase is taken as the resting phase. Criterion B: The cardiac cycle in which the rate of change of coronary artery cross-sectional area in adjacent time phases is less than a preset threshold is taken as the resting period; Criterion C: Use the stationary period specified in the scanning equipment protocol.

3. The calcification integral noise reduction method according to claim 1, characterized in that, Step 2 includes: Step 21: Read the ImagePositionPatient tag from the DICOM image header information of each time phase and extract the z-coordinate of each image layer; simultaneously read the spatial geometric information to confirm that the images of each time phase are in the same patient physical coordinate system; Step 22: For images at each time phase, calculate the common coverage area of ​​the z-coordinate range. The specific calculation formula includes: z_lo = max_i{min(Z_i)}; z_hi = min_i{max(Z_i)}; z_tol = 0.5 × SliceThickness; Where z_lo represents the lower bound of the common coverage area; Z_i represents the set of z coordinates of all slices of the image at the i-th time phase; z_hi represents the upper bound of the common coverage area; z_tol represents the half-slice thickness tolerance; and SliceThickness represents the slice thickness. For images at each time phase, select those that satisfy z_lo The physical mesh is constructed by extracting layers that are z_tol≤z≤z_hi+z_tol in order of z-coordinate. Step 23: If there are differences in the orientation matrix, pixel spacing, matrix size, or field of view of the images at different time phases, then the physical grid of the image at the anchor point time phase is used as the target grid to establish a mapping relationship from the reference time phase to the target grid. Step 24: When the common coverage area is less than a preset range threshold, individual images are missing layers, or there are inconsistent layer counts, execute any of the following fallback strategies: Backtracking strategy 1: Based on the complete z-coordinate range of the image at the anchor point, images at other time points are matched according to the closest physical plane; Backoff strategy 2: Take the minimum number of common layers with overlapping regions in all time phase images; Rollback strategy 3: For layers with missing edges, remove them only during the final integration; Rollback Strategy 4: When the image of any reference time phase cannot be mapped to the physical grid of the anchor time phase, the image of that reference time phase is discarded.

4. The calcification integral noise reduction method according to claim 1, characterized in that, The high-density voxels refer to voxels with a density value not lower than a preset density threshold T.

5. The calcification integral noise reduction method according to claim 4, characterized in that, The preset density threshold T ranges from 130 to 200 HU.

6. The calcification integral noise reduction method according to claim 1, characterized in that, Step 5 specifically includes: Step 51: Perform three-dimensional morphological dilation only on the candidate calcification mask of the reference phase; Step 52: Determine the displacement tolerance D_tol based on the maximum possible physiological displacement D_max between the image of the reference phase_i and the image of the anchor phase; Step 53: Calculate the number of discrete expansion iterations N_tol. The specific formula includes: N_tol=max(1,round(D_tol / PixelSpacing)); PixelSpacing refers to the axial pixel spacing.

7. The calcification integral noise reduction method according to claim 1, characterized in that, Step 7 specifically includes: Step 71: Suppose the Lth lesion is composed of a two-dimensional 8-connected component C_L from Mask_final, calculate the lesion area; the specific calculation formula is as follows: Area_L=PixelSpacing_x·PixelSpacing_y·#C_L; Where Area_L represents the lesion area; PixelSpacing_x represents the projection of the axial pixel spacing onto the x-axis; PixelSpacing_y represents the projection of the axial pixel spacing onto the y-axis; and #C_L represents the number of pixels contained in C_L. Step 72: Let h_L be the peak value of the image pixel HU in the time phase of the anchor point within C_L. The specific calculation formula is: h_L=max_{(x,y)∈C_L}HU_anchor(x,y); The Agatston density weights are calculated using the following formula: w(h)=1·1{130≤h<200}+2·1{200≤h<300}+3·1{300≤h<400}+4·1{h≥400}; Where w(h) represents the density weight value; 1{·} represents the indicator function; Step 73: Filter connected regions with lesion areas smaller than the preset minimum area A_min; combine with the coronary artery ROI, remove high-density connected regions located in non-coronary artery anatomical regions, and calculate the total modified Agatston integral, expressed by the following formula: S_total=Σ_{L∈ _valid}Score_L; Where S_total represents the total corrected Agatston integral value; _valid represents the set of lesions; Score_L represents the corrected Agatston integral value of the Lth lesion, calculated using the following formula: Score_L=Area_L·w(h_L)·F_thick; Where w(h_L) represents the density weight value of the Lth lesion; F_thick represents the thickness correction factor; Step 74: Construct five risk level functions based on the clinical Agatston risk grading guidelines.

8. The calcification integral noise reduction method according to claim 1, characterized in that, It also includes step 8: independently calculate the integral value on the images of each time phase, and calculate the coefficient of variation between time phases.

9. A coronary artery calcification integral noise reduction system, characterized in that, It includes a data acquisition module, a data processing module, and a result generation module; The data acquisition module is used to acquire images of the subject at multiple time phases during the quiescent period; the images are reconstructed images of PCCT ultrathin slices; one of the time phases is designated as the anchor time phase, and the remaining time phases are reference time phases; The data processing module includes a spatial positioning unit, an ROI generation unit, a mask generation unit, an expansion unit, a mask calculation unit, and an integral calculation unit; The spatial positioning unit is used to perform registration-free spatial positioning of images at each time phase based on the patient's physical coordinate system. The ROI generation unit is used to generate a coronary artery ROI on the image at the anchor point time phase, and to define the spatial range for candidate calcification extraction and final integral calculation of images at subsequent time phases. The mask generation unit is used to mark high-density voxels in images of each time phase within the coronary artery ROI, and generate candidate calcification masks for each time phase, including the candidate calcification mask Mask_anchor for the anchor point phase and the candidate calcification mask Mask_ref for the reference phase; the high-density voxels refer to voxels with a density value not lower than a preset density threshold T. The expansion unit is used to morphologically expand the candidate calcification mask of the reference time phase; The mask operation unit is used to perform a voxel-by-voxel OR operation on the candidate calcification masks of the expanded reference phase to obtain the confirmed mask. Then, the candidate calcification mask of the unexpanded anchor point phase is ANDed with the confirmed mask on a voxel-by-voxel basis to obtain the final calcification mask; The integral calculation unit performs calcification integral calculation on the image at the anchor point time phase based on the final calcification mask; The result generation module is used to send the processing results of the data processing module to external systems.

10. A computer-readable storage medium, characterized in that, After the program stored in the storage medium is executed by the processor, it can complete the steps of the calcification integral noise reduction method as described in any one of claims 1 to 8.

Images