A pathological wsi image discrete candidate screening method based on cycle consistency

Abstract

The application discloses a pathological WSI image discrete candidate screening method based on cycle consistency. In view of the problems that the existing method lacks closed-loop verification at the candidate level, the continuous domain cycle consistency is easy to fall into local minimum, and the weighted multiple graph lacks discrete multi-candidate competition, K candidates are generated in parallel by GPU and large-core CPU; a Markov random field (MRF) is constructed with image pairs as nodes; a pair potential energy ψ ab is constructed with a cycle consistency residual, and an optimal label is obtained by TRW-S; an abnormal node is triggered to be degraded and regenerated according to the sum S a of the pair potential energy; and a sparse adjacency matrix and an edge confidence vector are concatenated to a subsequent continuous domain global alignment module. The application is suitable for HE / IHC, fluorescence and Z-Stack modalities, and the alignment success rate is 99.2% and the TRW-S is only 4.2 seconds on an RK3588 platform with 5000 tiles.

Description

Technical Field

[0001] This invention relates to the field of digital pathological image processing and discrete combinatorial optimization technology, and in particular to a discrete candidate screening method for pathological WSI (Whole Slide Image) images based on constructing a Markov Random Field (MRF) with Cycle-Consistency residuals and solving it through Tree-Reweighted Message Passing (Sequential, TRW-S). This method runs on an embedded heterogeneous on-chip system composed of a Neural Processing Unit (NPU), a Graphics Processing Unit (GPU), a large-core CPU, and a small-core CPU. It serves as a preliminary discrete screening stage for the applicant's "A Global Alignment Method for Pathological WSI Images Based on Spectral SE(2) Synchronization and GNC Annealing" (hereinafter referred to as the "Main Patent"), which was filed on the same day. Background Technology

[0002] Digital pathology WSI scanners use objectives and linear or area array sensors to perform two-dimensional scanning of tissue sections, obtaining thousands to tens of thousands of local image tiles. These tiles are then stitched together according to their physical coordinate adjacency to form a complete full-field digital slide. However, due to sub-pixel positioning errors in mechanical scanning, uneven staining of tissue sections, and defocusing blurring at different focal planes, simple physical coordinate stitching produces significant misalignment, ghosting, and seams, failing to meet the dual requirements of millimeter-level field of view and micrometer-level resolution for pathological diagnosis.

[0003] Existing WSI stitching candidate generation and screening techniques are developing in three directions. The first type is the pairwise minimization plus ICP (Iterative Closest Point) refinement scheme, represented by TrakEM2 published by Preibisch et al. in Bioinformatics in 2009. It estimates an optimal relative transformation for each pair of overlapping tiles, including phase correlation, ECC (Enhanced Correlation Coefficient), SIFT (Scale-Invariant Feature Transform), and RANSAC (RANdom Sampling Consensus). Then, these pairwise transformations are adjusted using least squares on the global graph. Its limitations are as follows: when uneven staining, bubbles, folds, or focal plane differences exist in the overlapping area of ​​two tiles, pairwise estimation will yield an incorrect "confidence" transformation, which cannot be identified at the single-image level; the incorrect transformation is amortized to the entire image after least squares adjustment, causing global drift and local misalignment superposition; the incorrect transformation cannot be detected and removed subsequently unless manual review is introduced. Chinese patent CN110210447B (2020) discloses a pathological slide stitching method based on phase correlation, which can only handle pure translation models and is powerless for tile pairs with rotation or scaling, and also lacks a multi-candidate mechanism.

[0004] The second category is continuous optimization schemes that use cycle consistency as a regularization term. Represented by Christensen et al.'s "Image registration by minimizing a sorted vector of residuals," published in IEEE TMI in 2001, this approach assumes that the relative transformations of a set of images, when combined along closed loops, should equal a unit transformation, and smaller residuals indicate better overall consistency. Its limitations include: the potential energy is non-convex in the continuous domain, making it prone to getting trapped in local minima; when the initial candidate is incorrect, the cycle consistency constraint is unable to detect and replace it, only able to perform small-scale perturbations in the neighborhood of the incorrect candidate; and on large-scale WSI maps (thousands of image pairs, tens of thousands of closed loops), the computational cost of continuous domain joint optimization is extremely high, making it unacceptable for embedded devices.

[0005] The third category comprises joint optimization schemes based on weighted multigraphs, represented by US10,628,698 B1 (2020), US11,386,143 B2 (2022), and US12,167,161 B1 (2023) from Grundium, Finland. These schemes construct weighted multigraphs between image pairs, introduce virtual edges, and jointly minimize node transformations and edge weights using a linear constraint Σw = 1. Their limitations include: when systematic ambiguity exists in overlapping regions (such as periodic staining textures or tissue folds), pairwise estimations output multiple seemingly reasonable transformations, while the Grundium scheme only retains one "best" estimate for global optimization, leading to unresolved ambiguities that are amplified; a lack of a multi-candidate competition mechanism at the discrete level, with robustness entirely dependent on continuous domain weight softening; and the inability to directly perform label selection and cycle consistency joint verification on multiple candidates in the discrete domain.

[0006] Existing methods also have shortcomings in their implementation on embedded platforms. Solutions such as Christensen 2001 and Grundium rely on workstation-level CPUs and discrete GPUs, consuming more than 150 W and costing more than $5,000, making them unsuitable for integration into portable or desktop digital pathology scanners. The computational cost of continuous domain joint optimization increases exponentially with the number of closed loops, making it unscalable for WSIs of 5,000 to 20,000 tiles.

[0007] This invention addresses the aforementioned technical problems by providing a method for discrete candidate screening of pathological WSI images based on cycle consistency. It constructs an MRF using image pairs (not the original images) as nodes, constructs discrete pair potential energy using cycle consistency residuals, employs TRW-S for efficient solution on an embedded CPU, and combines anomaly node detection and candidate downgrading / regeneration mechanisms. On embedded platforms such as the Rockchip RK3588s (a Rockchip RK3588s embedded system-on-a-chip integrating 4 × ARM Cortex-A76 large cores, 4 × ARM Cortex-A55 small cores, an ARM Mali-G610 MP4 graphics processing unit, and 6 TOPS @INT8 neural network processing units), TRW-S convergence is achieved within 4.2 seconds for 5000 tile samples. The screening results are then cascaded to the continuous domain global alignment module of the main patent, resulting in an overall process time of less than 10 seconds. Summary of the Invention

[0008] To address the technical problems in the existing technologies, such as pairwise minimization without closed-loop verification, continuous domain cyclic consistency easily getting trapped in local minima and not being replaceable at the candidate level, weighted multigraphs lacking a discrete multi-candidate competition mechanism, and excessive power consumption and cost of embedded platforms, this invention provides a discrete candidate screening method for pathological WSI images based on cyclic consistency.

[0009] In a first aspect, the present invention provides a method for discrete candidate screening of pathological WSI images based on cycle consistency, comprising the following steps: S1: Collect N local image blocks (tiles) obtained from digital pathology full-view digital slide scanning, construct the original image G = (V, E) according to the physical coordinate adjacency relationship, and generate K candidate relative transformations T_ij(k) in parallel for each image pair (i, j) ∈ E, k = 1, 2, …, K, with a pairwise matching cost C_ij(k). S2: Construct a Markov random field (MRF), i.e., an image pair graph G_pair = (V_pair, E_pair), with image pairs as nodes and the relationship between image pairs sharing the same original image as edges. The label space of node a ∈ V_pair is {1, 2, …, K}, and the label l_a indicates that the image pair selects the l_a-th candidate transformation. S3: For each node a and each candidate label l_a, construct a univariate potential energy U_a(l_a) = C_ij(l_a)− min_k C_ij(k) + ε, where ε is a small positive number to avoid all zeros; S4: For each edge (a, b) in the image pair graph E_pair, construct the pair potential ψ_ab(l_a, l_b) = λ · ‖T_ij(l_a) · T_jk(l_b) − T_ik*‖_F² based on the cyclic consistency residual of the two image pairs sharing the same original image, where T_ik* is the expected composite transformation of the image pair (i, k), and λ is the pair potential weight coefficient; S5: Based on the univariate potential energy U_a(l_a) and the relative potential energy ψ_ab(l_a, l_b), define the total energy of the MRF as E(L) = Σ_a U_a(l_a) + Σ_(a,b) ψ_ab(l_a, l_b), where L is the label assigned to all nodes; S6: The sequence tree reweighted message passing algorithm TRW-S is used to solve the discrete energy minimization problem of the total energy E(L) of the MRF, and the optimal candidate label l_a* and the corresponding candidate transformation T_ij(l_a*) for each image pair are output. S7: For each node a after TRW-S convergence, calculate the sum of potential energy under its optimal label, S_a = Σ_b ψ_ab(l_a*, l_b*). When S_a is greater than the threshold θ_high, mark the node as an abnormal node and perform candidate downgrade regeneration and local iteration on the node and its 2-hop neighborhood. When S_a is still greater than θ_high after local iteration, delete the image pair from the original image G. S8: Output candidate filtering results, which include a filtered candidate index table, a sparse adjacency matrix à and an edge confidence vector w, where w_a = exp(−S_a / σ²), which serves as the input to the subsequent continuous domain global alignment module; S9: The global task scheduler of the embedded on-chip system delivers the sparse adjacency matrix Ã, candidate transformation T_ij(l_a*) and edge confidence vector w to the subsequent continuous domain global alignment module through shared memory, and performs subsequent global alignment in a cascaded manner; S10: The global task scheduler performs heterogeneous resource scheduling on S1 to S9, allocating tasks with different computational characteristics to the neural network processing unit (NPU), graphics processing unit (GPU), large-core CPU, and small-core CPU for execution.

[0010] Furthermore, the parallel generation of K candidate relative transformations T_ij(k) for each image pair (i, j) ∈ E in S1 includes: the GPU performing phase correlation in parallel FFT (Fast Fourier Transform) to obtain phase correlation candidates; the CPU performing SIFT and RANSAC to obtain K-2 geometric verification candidates; the CPU performing ECC mutual information estimation to obtain mutual information candidates; merging and deduplicating to retain K candidates; the value of K ranges from 3 to 8.

[0011] Furthermore, the construction of the potential energy ψ_ab(l_a, l_b) described in S4 includes: for forward sharing, ψ_ab = λ · ‖T_ij(l_a) · T_jk(l_b) − T_ik*‖_F²; for reverse sharing, ψ_ab = λ · ‖T_ij(l_a) ·T_ik(l_b)⁻¹ − T_jk*‖_F²; T_ik* takes the candidate with the minimum cost when (i, k) ∈ E, and takes the unit transformation I when (i,k) ∉ E; λ takes values ​​from 0.5 to 2.0.

[0012] Furthermore, the TRW-S discrete energy minimization solution described in S6 includes: decomposing the MRF graph into multiple subtrees; dividing node traversal into blocks based on image physical coordinates, with each block running in parallel using OpenMP (Open Multi-Processing) multithreading and between blocks in a topological order sequence; vectorizing the K × K potential energy matrix of each edge using SIMD (ARM Single Instruction Multiple Data Stream instruction set, used for vectorized message update core loop) instructions; and stopping iteration when the relative energy change between two adjacent iterations is less than the convergence threshold τ or the maximum number of iterations iter_max is reached.

[0013] Furthermore, the abnormal node detection and candidate downgrading and regeneration described in S7 includes: calculating the sum of potential energy S_a for each node; marking an abnormal node when S_a is greater than the threshold θ_high; downgrading and regenerating the candidate set of abnormal nodes (expanding the SIFT matching radius, reducing the RANSAC threshold, and adding LOG (Logarithm) polar coordinate candidates); re-executing the local iterations of S3 to S6 only for the 2-hop neighborhood of the abnormal node; and deleting the image pair from the original image G when S_a is still greater than θ_high after the local iterations.

[0014] Furthermore, the edge confidence vector w_a = exp(−S_a / σ²) mentioned in S8, where σ ranges from 1.0 to 3.0; w_a serves as the initial robust weight for GNC annealing in the subsequent continuous domain global alignment module.

[0015] Furthermore, the cascaded execution of subsequent global alignment described in S9 includes: zero-copy delivery of data via DMA-BUF (Linux kernel DMA buffer sharing framework for zero-copy data transfer across processing units); subsequent continuous domain global alignment module using à and T_ij(l_a*) as input to execute spectrum SE(2) to synchronously calculate the initial global pose; using w as initial weight to perform GNC annealing to remove residual outliers; outputting the final global pose and entering multi-band fusion.

[0016] Furthermore, the heterogeneous resource scheduling described in S10 includes: the small-core CPU performing tile decoding, physical coordinate management, and global task scheduling; the GPU performing FFT phase correlation using OpenCL (Open Computing Language); the large-core CPU performing MRF construction, TRW-S message passing, and abnormal node detection using OpenMP; data transfer between processing units via DMA-BUF zero-copy; and dynamic frequency adjustment based on load using DVFS (Dynamic Voltage and Frequency Scaling).

[0017] Furthermore, the method is applied to three modalities: bright-field HE (Hematoxylin & Eosin) staining, immunohistochemical IHC (Immunohistochemistry) staining, and fluorescence multi-channel and image stack Z-Stack three-dimensional tomography. The only difference between the different modalities is the replacement of the local feature sources generated by the candidate in S1.

[0018] Secondly, the present invention provides a discrete candidate screening device for pathological WSI images based on cycle consistency, comprising: a candidate generation module, an image pair graph construction module, a potential energy construction module, a discrete energy minimization module, an abnormal node detection and degradation module, a candidate screening result output module, a cascaded execution module, and a heterogeneous scheduling module.

[0019] Thirdly, the present invention provides an electronic device, including a processor and a memory, wherein the memory stores a computer program, and the processor executes the program to implement the steps of the above-described method.

[0020] Fourthly, the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the steps of the above-described method.

[0021] Compared with the prior art, the present invention has the following beneficial effects: (1) By constructing an MRF with image pairs (instead of the original image) as nodes, the label of each node is the candidate transformation index of the image pair, so that the "candidate screening" problem can be expressed as a discrete energy minimization problem, which is fundamentally different from the continuous domain schemes such as Grundium US10,628,698 B1 with the original image as nodes. Incorrect candidates can be directly eliminated at the global closed-loop level rather than being disturbed in the neighborhood. (2) By constructing a pair of potential energy ψ_ab(l_a, l_b) = λ · ‖T_ij(l_a)· T_jk(l_b) − T_ik*‖_F² based on the cycle consistency residual, the cycle consistency is transformed from the continuous regularization term of Christensen 2001 to the discrete MRF pair of potential energy for the first time, so that the elimination of erroneous candidates is hard-selected at the candidate level, avoiding the local minima of continuous domain non-convex optimization; (3) By using the TRW-S algorithm to minimize MRF energy, the sequence tree reweighted message passing converged 5000 tile samples in 4.2 seconds on the RK3588s embedded CPU. This is the first engineering deployment of TRW-S in the WSI splicing field. (4) By dividing nodes into blocks according to the physical coordinates of the image, adding OpenMP multi-threaded parallelism and NEON SIMD vectorized message passing, the monotonic convergence of TRW-S and the utilization of embedded multi-core are taken into account. The message update core loop is accelerated by 3.2 times in actual tests. (5) Through the abnormal node detection and candidate downgrade regeneration mechanism, the abnormal node is automatically identified based on the criterion that the sum of potential energy S_a is greater than the threshold θ_high. The candidate downgrade regeneration and local iteration are triggered in the 2-hop neighborhood of the node, so that the method is adaptive to typical failure scenarios of pathological sections such as tissue folding, bubbles, and focal plane drift, without the need for manual intervention. (6) By constructing a pair-graph, candidate screening can be extended to tens of thousands of tiles for the first time. The complexity of TRW-S is O(P · K² · iter) ≈ O(N · K² · iter), which is scalable for 5,000 to 20,000 tiles. Compared with the continuous domain schemes such as Christensen 2001, the complexity increases exponentially with the number of closed loops. (7) By concatenating the screening results (sparse adjacency matrix Ã, candidate transformation T_ij(l_a*), edge confidence vector w) to the continuous domain global alignment module of the main patent, the convergence domain of the main patent is significantly expanded and its convergence is accelerated. The overall time for the entire process of 5000 tiles on RK3588s is less than 10 seconds, and the alignment success rate of pathological sections with "uneven staining and periodic texture" is increased from 87% of the Grundium-type scheme to 99.2%. Attached Figure Description

[0022] Figure 1 is an overall flowchart of the pathological WSI image discrete candidate screening method based on cycle consistency in an embodiment of the present invention, showing ten steps from S1 multi-source candidate generation to S10 heterogeneous resource scheduling, and S7 abnormal node detection. a The relationship between the two comparisons with the threshold θ_high and the branch processing such as downgrade regeneration and edge deletion of the original image.

[0023] Figure 2 is a schematic diagram of the parallel scheduling of multi-source candidate generation in S1 in this embodiment of the invention. It shows that the input tilesi and j are sent to the GPU and the CPU for parallel processing via libtiff image decoding and DMA-BUF zero-copy by the small-core CPU Cortex-A55: the GPU Mali-G610 performs two-dimensional FFT, cross-power spectrum phase addition, IFFT peaking and log-polar rotation estimation to output phase-correlated candidate T_ij⁽¹⁾; the large-core CPU Cortex-A76 obtains K-2 geometric candidates from SIFT keypoints through RANSAC outlier removal, and obtains mutual information candidate T_ij⁽ by ECC mutual information grid scanning. 5 The three candidate classes are sorted after deduplication and merged with a translation difference of less than 0.5 pixels to obtain K candidate T_ij⁽ᵏ⁾ (k = 1…K) in parallel pipeline relationship.

[0024] Figure 3 is a schematic diagram of the potential energy construction of the image pair graph MRF and the cyclic consistency pair in embodiments of the present invention (S2 to S5). It shows that the adjacency relationship of tiles i, j, k in the original graph G = (V, E) is mapped to the nodes a = (i, j), b = (j, k), c = (i, k) in the image pair graph G_pair = (V_pair, E_pair), with the node label space {1, 2, …, K}, and the two pairs of potential energy ψ: forward sharing (sharing the middle image j) and backward sharing (sharing the left image i). a ᵦ Construction situation (forward ψ) a ᵦ = λ · ‖T_ij(l a ) · T_jk(l_b) − T_ik*‖²;Inverse ψ a ᵦ = λ · ‖T_ij(l a ) · T_ik(l_b)⁻¹ − T_jk*‖²), the value rules for the expected composite transformation T_ik* ((i, k) ∈ E take argmin C_ik, (i, k) ∉ E take the unit transformation I), univariate potential energy U a (l a ) = C_ij(la ) − min k C_ij(k) + ε and the potential energy ψ a ᵦ The total energy E(L) is equal to Σ a U a (l a ) + Σ₍ a ,ᵦ₎ ψ a ᵦ(l a The process of , l_b).

[0025] Figure 4 is a schematic diagram of S6 TRW-S solving plus S7 abnormal node detection plus degradation and regeneration in an embodiment of the present invention. It shows the core process of the TRW-S solver, including tree decomposition plus OpenMP block parallelism, NEON SIMD message update, and monotonically decreasing energy to convergence, outputting the optimal label l_a* and candidate transformation T_ij(l_a*); and the sum of potential energy S in S7. a = Σᵦ ψ a The process involves comparing ᵦ(l_a*, l_b*) with the threshold θ_high, downgrading and regenerating abnormal node candidates, iterating through the local neighborhood of 2 hops, and then re-determining whether to restore the node to a normal node or delete the edge from the original graph G.

[0026] Figure 5 is a resource scheduling architecture diagram of the embedded heterogeneous on-chip system in an embodiment of the present invention. It shows the global task scheduler and DVFS dynamic frequency adjustment mechanism running on the small-core CPU, as well as the data path of the three types of processing units, namely the small-core CPU Cortex-A55 (image decoding, physical coordinate management, global task scheduling, file I / O), GPU Mali-G610 (OpenCL FFT phase correlation, parallel candidate generation stage, data path to candidate merging) and large-core CPU Cortex-A76 (MRF construction, TRW-S message passing, abnormal node detection, SIFT plus RANSAC plus ECC), which are jointly executed through the Linux DMA-BUF zero-copy shared memory layer and L1 cache LUT 1024 items.

[0027] Figure 6 is a schematic diagram of the cascaded execution of the discrete screening stage before the main patent in the embodiment of the present invention. From left to right, it shows the output of the discrete candidate screening stage of the present invention: sparse adjacency matrix Ã, candidate transformation T_ij(l_a*) and edge confidence vector w. These are delivered to the global alignment module of the continuous domain of the main patent via Linux DMA-BUF zero-copy shared memory. The module then sequentially executes the spectrum SE(2) to synchronously calculate the initial global pose, GNC annealing with Geman-McClure robust kernel, remove residual outliers, finally solve the global pose, and perform multi-band fusion. Finally, it outputs the pathological WSI full-view digital slice. Detailed Implementation

[0028] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features described in these embodiments can be combined with each other.

[0029] Example 1: Bright-field HE staining WSI discrete candidate screening (5000 tiles, no abnormal nodes) Referring to Figure 1, this embodiment is implemented on a Rockchip RK3588s embedded system-on-a-chip (4 Cortex-A76 large cores @ 2.4 GHz + 4 Cortex-A55 small cores @ 1.8 GHz + Mali-G610 MP4 GPU + 6 TOPS NPU, power consumption less than 30 W), running an embedded Linux 6.1 kernel and a user-space stitching program. The sample under test is a liver tissue HE-stained section, which is scanned by a 40x objective lens and a 4096 × 1 linear array sensor to obtain N = 5000 tiles, each tile being 4096 × 4096 pixels in size. The physical coordinates are given by the encoder of the scanning station. The adjacency relationship is constructed as the original image G = (V, E), |V| = 5000, |E| = 9897 (including four adjacent and a small number of diagonal adjacent).

[0030] S1 (Generating the set of candidate transforms for image pairs): For each image pair (i, j) ∈ E, generate K = 5 candidate relative transforms T_ij(k) in parallel, k = 1, 2, …, 5, along with a pairwise matching cost C_ij(k). See Figure 2. This step is executed by the candidate generation module under the coordination of the heterogeneous scheduling module. S1.1 uses a small-core CPU (Cortex-A55) with a libtiff thread to decode the TIFF and physical coordinate metadata of the tile in parallel, and sends the YUV420 frame to the dedicated memory heap of the GPU and the large-core CPU in a zero-copy manner through the Linux kernel DMA-BUF framework. S1.2 Phase correlation is performed by GPU (Mali-G610 MP4) in parallel using OpenCL: 2D FFT is performed on the grayscale images of the two tiles, the cross power spectrum phase is taken, and then IFFT is performed. The peak position is the translation component of the phase correlation candidate T_ij(1); rotation and scaling are obtained by phase correlation after log-polar transformation and merged into T_ij(1). The cost C_ij(1) = -peak height; S1.3 The large-core CPU (Cortex-A76) performs SIFT plus RANSAC in 4 threads using OpenMP to obtain K-2 = 3 geometric verification candidates T_ij(2), T_ij(3), T_ij(4), with a cost C_ij(2..4) = the reciprocal of the number of interior points + 0.01 · residual variance; S1.4 The enhanced correlation coefficient ECC mutual information estimation of the grayscale histograms of the two tiles is calculated by the big core CPU. A 9 × 9 grid scan is performed on the translation parameter in the ±8 pixel neighborhood. The grid point with the largest ECC is used as the mutual information candidate T_ij(5), and the cost C_ij(5) = -ECC. S1.5 Merge K = 5 candidates and remove duplicates (translation component differences less than 0.5 pixels and rotation differences less than 0.2° are considered as the same candidate, and the one with the lower cost is retained), and output the candidate set {T_ij(1), …, T_ij(5)} and the cost vector {C_ij(1), …, C_ij(5)}.

[0031] In the formula, T_ij(k) ∈ SE(2) is a 3 × 3 homogeneous transformation matrix that maps image j to the coordinate system of image i, with 3 degrees of freedom (translation t_x, t_y and rotation θ); C_ij(k) is a non-negative real number with a dimension of 1, and the smaller the value, the better the candidate quality; K = 5 is the number of candidates retained for each image pair.

[0032] S2 (Image Pair Graph MRF Construction): Referring to Figure 3, the image pair graph construction module constructs an image pair graph G_pair = (V_pair, E_pair) with image pairs (i, j) ∈ E as nodes, where |V_pair| = M = 9897. On the image pair graph, if two image pairs (i, j) and (j, k) share the same intermediate image j (i.e., the endpoint of the former is the starting point of the latter), then an edge (a, b) ∈ E_pair is connected between them; if two image pairs (i, j) and (i, k) share the same left image i (reverse sharing), an edge is also connected; |E_pair| = P ≈ 2 · M = 19794; the label space of node a ∈ V_pair is {1, 2, 3, 4, 5}, and the label l_a indicates that the image pair selects the l_ath candidate transformation.

[0033] In the formula, M is the total number of image pairs, which is equal to the number of edges in the original image |E|; P is the number of edges in the image pair graph, that is, the total number of "sharing relationships" between two image pairs that share the same original image, typically P ≈ 3M to 5M.

[0034] S3 (Unary Potential Energy Construction): The potential energy construction module constructs univariate potential energy for each node a (corresponding to the image pair (i, j)) and each candidate label l_a. U_a(l_a) = C_ij(l_a) − min_k C_ij(k) + ε; In the formula, ε is a small positive number to avoid all zeros, and ε = 0.01 is taken; C_ij(l_a) is taken from the pairwise matching cost of S1; U_a(l_a) has a dimension of 1 and its minimum value is equal to ε; the univariate potential energy reflects "the matching quality of the candidate in the image to itself", and using it alone cannot identify erroneous candidates that have "small pairwise costs but contradict the global closed loop".

[0035] S4 (Polypotential Energy Construction Based on Cyclic Consistency Residuals): Referring to Figure 3, the potential energy construction module constructs a pair of potential energy for each edge (a, b) in E_pair, in two cases: Scenario 1 (Forward Sharing): a = (i, j), b = (j, k) sharing intermediate image j, with potential energy of ψ_ab(l_a, l_b) = λ · ‖T_ij(l_a) · T_jk(l_b) − T_ik*‖_F²; Scenario 2 (Reverse Sharing): a = (i, j), b = (i, k) share the left graph i, with potential energy of ψ_ab(l_a, l_b) = λ · ‖T_ij(l_a) · T_ik(l_b)⁻¹ − T_jk*‖_F²; In the formula, ||·||_F is the Frobenius norm of the matrix, which for a 3 × 3 matrix M is ||M||_F² = Σ_pΣ_q M_pq²; T_ij(l_a) and T_jk(l_b) are candidate transformations of image pairs (i, j) and (j, k) under labels l_a and l_b, respectively, with 3 degrees of freedom; T_ik* is the expected composite transformation of image pair (i, k). When (i, k) ∈ E, i.e., i and k are adjacent image pairs, the candidate T_ik(argmin_l C_ik(l)) with the minimum pairwise matching cost in the candidate set of (i, k) is selected. When (i, k) ∉ E, the unit transformation I is selected; λ is the potential weight coefficient, with a value range of 0.5 to 2.0. In this embodiment, λ = 1.0 is taken; ψ_ab(l_a, l_b) The dimension is 1, and the smaller the value, the better the cycle consistency of the label combination (l_a, l_b).

[0036] The physical meaning of the potential energy is as follows: if the composite of two images sharing the same original image does not match the "expected composite" after combining their respective selected candidates, a high penalty score is given; otherwise, a zero penalty score is given. This is the core difference between this invention and the prior art—transforming the cycle consistency from the continuous regularization term of the Christensen 2001 formula to the discrete MRF pair potential energy, so that erroneous candidates are directly eliminated at the global closed-loop level.

[0037] S5 (Definition of Total Energy in MRF): The total energy of the entire MRF is... E(L) = Σ_a U_a(l_a) + Σ_(a,b) ∈ E_pair ψ_ab(l_a, l_b); In the formula, L = {l_a | a ∈ V_pair} is the label assignment for all nodes; E(L) has a dimension of 1, and the smaller the value, the better the label assignment matches the global loop consistency; the goal of this invention is to solve L* = argmin_L E(L).

[0038] S6 (TRW-S Discrete Energy Minimization Solution): Referring to Figure 4, the discrete energy minimization module uses the TRW-S algorithm proposed by Kolmogorov PAMI 2006 to minimize E(L): S6.1 Decompose the MRF graph into multiple subtrees according to the maximum spanning tree, with each edge belonging to at least one subtree; S6.2 Node traversal is divided into blocks based on the physical coordinates of the image (each block is 16 × 16 nodes). Within each block, OpenMP multi-threaded parallelism is used (4 threads run on the Cortex-A76 core), and between blocks, the topological order is followed to ensure monotonically decreasing energy. S6.3 stores the 5 × 5 potential energy matrices for each edge in row-major order and uses ARM NEON single instruction multiple data stream SIMD instructions (vld1q_f32 × 4 loads 16 float32s at a time) to vectorize the message update core loop; S6.4 When the relative energy change between two adjacent iterations (E_(t-1) − E_t) / E_(t-1) is less than the convergence threshold τ or the maximum number of iterations iter_max is reached, the iteration stops and the optimal label l_a* and the corresponding candidate transformation T_ij(l_a*) of each node are output; in this embodiment, τ = 1e-3 and iter_max = 50 are taken.

[0039] In the formula, E_t is the total energy of MRF after the t-th iteration; τ is the relative convergence threshold, ranging from 1e-4 to 1e-3; iter_max is the maximum number of iterations, ranging from 30 to 100; in this embodiment, 5000 tile samples converged within 30 iterations, and the TRW-S stage took 4.2 seconds.

[0040] S7 (Abnormal Node Detection and Candidate Degradation Regeneration): Referring to Figure 4, after TRW-S converges, the abnormal node detection and degradation module checks the sum of the potential energies of the final label l_a* for each node a: S_a = Σ_(b ∈ N(a)) ψ_ab(l_a*, l_b*); In the formula, N(a) is the set of neighbors of node a in the image pair graph; S_a is dimensionless and the equivalent displacement is in square pixels.

[0041] S7.1 If S_a is less than or equal to the threshold θ_high (θ_high = 5.0 in this embodiment), then the node is considered a normal node and no degradation is triggered; S7.2 If S_a is greater than θ_high, the node is marked as an "abnormal node," meaning that no matter which candidate is selected, the cyclic consistency of the neighborhood cannot be satisfied, and the overall candidate quality of the image pair is poor. The processing strategy is to "downgrade and regenerate" the candidate set of abnormal nodes (increase the SIFT matching radius by 1.5 times, reduce the RANSAC threshold, and add 2 LOG polar coordinate candidates). S7.3 only re-executes the local iterations of S3 to S6 for the 2-hop neighborhood of the abnormal node (i.e., the node and its direct neighbors and the neighbors of its neighbors) to avoid recalculating the entire graph; S7.4 If S_a is still greater than θ_high after local iteration, then the candidate of the image pair is discarded, the edge is deleted from the original image G, and the corresponding position is set to zero in the sparse adjacency matrix Ã.

[0042] In the formula, θ_high is the threshold for judging abnormal nodes, with a value range of 1.0 to 5.0, and the unit is the equivalent displacement of square pixels; in this embodiment, no abnormal nodes were triggered in the 5000 tile bright field HE sample (all nodes S_a are less than θ_high).

[0043] S8 (Candidate Screening Result Output): The candidate screening result output module outputs three documents to the contiguous field global alignment module of the main patent. (1) Candidate index table after filtering: the optimal label l_a* and the corresponding candidate transformation T_ij(l_a*) for each image pair; (2) Sparse adjacency matrix Ã: retain image pairs that are not empty (after deleting abnormal nodes) for the main patent's spectrum SE(2) synchronization; (3) Edge confidence vector w: w_a = exp(−S_a / σ²), in this embodiment σ = 2.0; w_a has a dimension of 1 and a value range of 0 to 1. The closer w_a is to 0, the less reliable the candidate of the image pair is, and the closer w_a is to 1, the more reliable it is; w is used as the initial robust weight of GNC annealing in the main patent.

[0044] In the formula, σ is a temperature parameter, ranging from 1.0 to 3.0, which controls the sensitivity of the confidence vector w to the sum of potential energies S_a.

[0045] S9 (Cascaded execution with the global alignment module of the main patent's continuous domain): See Figure 6. This method serves as a pre-discretionary screening stage for the main patent and is executed in cascades: S9.1 After this method is completed (S1 to S8), the three data sets Ã, T_ij(l_a*), and w are delivered to the contiguous domain global alignment module of the main patent through the zero-copy shared memory of the Linux kernel DMA-BUF framework; S9.2 The main patent uses à and T_ij(l_a*) as input to execute the spectrum SE(2) synchronously to obtain the global pose initial value {(R_i, t_i)}; S9.3 The main patent performs GNC annealing with w as the initial weight, further eliminating residual outliers in the continuous domain; S9.4 outputs the final global pose from the main patent and enters multi-band fusion.

[0046] After cascading, the entire process of aligning 5000 tiles on the RK3588s takes approximately 9.5 seconds (4.2 seconds for this method 2B + 5.3 seconds for the main patent 2A), which meets the engineering goal of "completing global alignment of 5000 tiles in less than 10 seconds".

[0047] S10 (Embedded Heterogeneous Resource Scheduling): As shown in Figure 5, the global task scheduler of the heterogeneous scheduling module (running on a small-core CPU) performs heterogeneous resource scheduling for S1 to S9. (1) The small-core CPU (Cortex-A55 × 4) performs tile image decoding, physical coordinate management and global task scheduling; (2) The GPU (Mali-G610 MP4) performs parallel computation of the FFT phase correlation and candidate generation stage in OpenCL (S1.2). (3) The large-core CPU (Cortex-A76 × 4) executes MRF construction, TRW-S message passing and abnormal node detection (S2 to S7) in OpenMP multi-threaded mode. (4) Tensor, image and matrix data transfer between processing units is performed in a zero-copy manner through the Linux kernel DMA-BUF framework; (5) The voltage and frequency of each processing unit are dynamically adjusted by the dynamic voltage and frequency adjustment DVFS mechanism according to the current stage load: GPU is fully loaded and CPU frequency is reduced during the candidate generation stage; CPU is fully loaded and GPU frequency is reduced during the TRW-S stage; zero-copy memory is used in the output stage to avoid repeated copying; when the abnormal node is downgraded and regenerated (S7) is triggered, only the affected A76 core is woken up, and the remaining cores maintain low power consumption.

[0048] Total process time: Candidate generation 1.8 s + MRF construction and potential energy calculation 0.9 s + TRW-S solution 1.5 s + abnormal node detection 0.3 s + output 0.2 s + other 0.4 s = 5.1 s (excluding the main patent 2A part). After adding the main patent 2A of about 4.4 s, the overall time is 9.5 s, which meets the target.

[0049] Example 2: Bright-field HE staining WSI discrete candidate screening (including downgrading and regenerating anomalous nodes) The sample and hardware platform in this embodiment are the same as those in Embodiment 1. The difference is that there are local oil droplets and bubbles in the objective lens and tissue folds during the scanning process, which causes S7 to detect 320 abnormal nodes in M ​​= 9897 image pairs (S_a> θ_high = 5.0).

[0050] S7 (Abnormal Node Detection and Candidate Degradation Regeneration): See Figure 4, executed by the abnormal node detection and degradation module: S7.1 For each node a, calculate the sum of potential energy under its optimal label: S_a = Σ_(b ∈ N(a))ψ_ab(l_a*, l_b*); If the S_a of 320 image pairs exceeds the threshold θ_high = 5.0 (maximum 18.7), they are marked as abnormal nodes. S7.2 Perform downgraded regeneration on the candidate set of the 320 anomalous nodes: increase the SIFT matching radius from ±8 pixels to ±12 pixels, increase the RANSAC threshold from 1.5 pixels to 2.5 pixels, and add 2 LOG polar coordinate candidates; S7.3 only re-executes the local iterations of S3 to S6 for the 2-hop neighborhood (a total of 1847 image pairs) of the 320 abnormal nodes to avoid recalculating the entire image; After local iteration in S7.4, the S_a of 287 image pairs drops below θ_high, restoring them to normal node state; the remaining 33 image pairs still have S_a greater than θ_high, so they are discarded as a whole, deleted from the original image G, and set to zero in the sparse adjacency matrix Ã.

[0051] In the formula, θ_high is the threshold for judging abnormal nodes; in this embodiment, 320 / 9897 ≈ 3.2% of the image pairs in the out-of-sample were triggered to be downgraded, and finally 33 / 9897 ≈ 0.33% of the image pairs were discarded. The main patent 2A can still complete global alignment on the remaining 9864 image pairs, with a final alignment success rate of 99.2% (a significant improvement compared to the 87% of the Grundium-type scheme).

[0052] The local iteration took 0.8 seconds. After adding the 5.1 seconds of Example 1, the total time of this method 2B was 5.9 seconds. After adding the main patent 2A for about 4.4 seconds, the overall time was 10.3 seconds (slightly exceeding the 10-second threshold by 0.3 seconds, but still 29 times faster than the workstation solution of 5 minutes).

[0053] Example 3: Fluorescence Multichannel and Z-Stack 3D Tomography Mode Extension Referring to Figure 6, this embodiment illustrates the differences when the method of the present invention is applied to three modalities: bright-field HE and IHC, fluorescence multichannel and Z-Stack three-dimensional tomography. The three modalities share the discrete candidate screening backends S2 to S9, only replacing the local feature sources of S1: (1) Bright field HE and IHC staining modalities: The light intensity changes caused by staining are gradual and the tissue texture is rich. SIFT and normalized cross-correlation NCC are used as local feature sources. When generating candidates, K = 5 candidates are generated from GPU FFT phase correlation, large-core CPU SIFT plus RANSAC (3 candidates), and large-core CPU ECC mutual information (1 candidate). (2) Fluorescence multi-channel mode: There is cross-color and background autofluorescence between channels, but the signal within each channel is sparse and there are few SIFT feature points. SIFT and phase correlation combination is used as local feature source; when generating candidates, K = 5 candidates are generated from GPU FFT phase correlation (including log-polar rotation estimation), large-core CPU SIFT plus RANSAC (2), and large-core CPU ECC mutual information (2). (3) Z-Stack 3D layer scanning mode: M z layers are acquired at each xy position, and the focal plane of each layer is different; SIFT and Sparse Microscopy Localization (SML) focal measurement are used as local feature sources; when generating candidates, K = 5 candidates are generated from large-core CPU SIFT plus RANSAC (3 candidates) and large-core CPU SML focal measurement (2 candidates, the layer with the highest focal peak is used as the representative frame to participate in the stitching). The three modalities share the discrete candidate screening backends S2 to S9, demonstrating the modal scalability of the method described in this invention. Device Examples

[0054] Referring to Figures 1 to 6, a second aspect of the present invention provides a discrete candidate screening device for pathological WSI images based on cycle consistency, comprising: The candidate generation module is used to acquire N tiles and generate K candidate relative transformations for each image pair; The image-pair graph construction module is used to construct Markov random fields (MRFs) using image pairs as nodes. The potential energy construction module is used to construct unary potential energy based on pairwise matching cost and pair potential energy based on cycle consistency residual; The discrete energy minimization module is used to solve for the optimal candidate label for the total energy of the MRF using the TRW-S algorithm; The abnormal node detection and degradation module is used to detect abnormal nodes and perform candidate degradation regeneration and local iteration on the node and its 2-hop neighborhood; The candidate selection result output module is used to output the sparse adjacency matrix Ã, the candidate transformation T_ij(l_a*), and the edge confidence vector w; A cascaded execution module is used to deliver candidate screening results to the contiguous domain global alignment module of the main patent via shared memory. The heterogeneous scheduling module is used to allocate the tasks of the various modules to the neural network processing unit (NPU), graphics processing unit (GPU), large-core CPU, and small-core CPU for execution.

[0055] The specific implementation methods of each module correspond one-to-one with S1 to S10 of the method described in the first aspect, and will not be repeated here. Equipment Examples

[0056] A third aspect of the present invention provides an electronic device, including a processor and a memory, wherein the memory stores a computer program, and the processor executes the program to implement steps S1 to S10 of the method described in the first aspect. In practical applications, the electronic device can be implemented using embedded heterogeneous system-on-a-chip such as Rockchip RK3588s, Horizon Robotics Sunrise 5, and NVIDIA Jetson Orin Nano, and is not limited to any specific brand or model. Medium Examples

[0057] A fourth aspect of the present invention provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements steps S1 to S10 of the method described in the first aspect. The computer-readable storage medium includes, but is not limited to, ROM, RAM, eMMC, UFS, SD card, SSD, HDD, CD-ROM, DVD, optical storage, magnetic storage, etc.

[0058] The above description is merely an embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the present invention by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principle of the present invention should be included within the scope of the claims of the present invention.

Claims

1. A method for discrete candidate screening of pathological WSI images based on cycle consistency, characterized in that, Includes the following steps: S1: Collect N local image blocks (tiles) obtained from digital pathology full-view digital slide scanning, construct the original image G = (V, E) according to the physical coordinate adjacency relationship, and generate K candidate relative transformations T_ij(k) in parallel for each image pair (i, j) ∈ E, k = 1, 2, …, K, with a pairwise matching cost C_ij(k). S2: Construct a Markov Random Field (MRF), i.e., an image pair graph G_pair = (V_pair, E_pair), with image pairs as nodes and the relationship between image pairs sharing the same original image as edges. The label space of node a ∈ V_pair is {1, 2, ..., K}, and the label l_a indicates that the image pair selects the l_a-th candidate transformation. S3: For each node a and each candidate label l_a, construct a univariate potential energy U_a(l_a) = C_ij(l_a) −min_k C_ij(k) + ε, where ε is a small positive number to avoid all zeros; S4: For each edge (a, b) in the image pair graph E_pair, construct the pair potential ψ_ab(l_a, l_b) = λ · ‖T_ij(l_a) · T_jk(l_b) − T_ik*‖_F² based on the cyclic consistency residual of the two image pairs sharing the same original image, where T_ik* is the expected composite transformation of the image pair (i, k), and λ is the pair potential weight coefficient; S5: Based on the univariate potential energy U_a(l_a) and the relative potential energy ψ_ab(l_a, l_b), define the total energy of the MRF as E(L) = Σ_a U_a(l_a) + Σ_(a,b) ψ_ab(l_a, l_b), where L is the label assigned to all nodes; S6: The sequence tree reweighted message passing algorithm TRW-S is used to solve the discrete energy minimization problem of the total energy E(L) of the MRF, and the optimal candidate label l_a* and the corresponding candidate transformation T_ij(l_a*) for each image pair are output. S7: For each node a after TRW-S convergence, calculate the sum of potential energy under its optimal label, S_a = Σ_b ψ_ab(l_a*, l_b*). When S_a is greater than the threshold θ_high, mark the node as an abnormal node and perform candidate downgrade regeneration and local iteration on the node and its 2-hop neighborhood. When S_a is still greater than θ_high after local iteration, delete the image pair from the original image G. S8: Output candidate filtering results, which include a filtered candidate index table, a sparse adjacency matrix à and an edge confidence vector w, where w_a = exp(−S_a / σ²), which serves as the input to the subsequent continuous domain global alignment module; S9: The global task scheduler of the embedded on-chip system delivers the sparse adjacency matrix Ã, candidate transformation T_ij(l_a*) and edge confidence vector w to the subsequent continuous domain global alignment module through shared memory, and performs subsequent global alignment in a cascaded manner; S10: The global task scheduler performs heterogeneous resource scheduling on S1 to S9, allocating tasks with different computational characteristics to the neural network processing unit (NPU), graphics processing unit (GPU), large-core CPU, and small-core CPU for execution.

2. The method according to claim 1, characterized in that, The parallel generation of K candidate relative transformations T_ij(k) for each image pair (i, j) ∈ E as described in S1 further includes: The GPU performs phase correlation using a parallel Fast Fourier Transform (FFT) to obtain phase correlation candidates. The large-core CPU performs Scale Invariant Feature Transform (SIFT) and Random Sample Consensus (RANSAC) to obtain K-2 geometric verification candidates. The large-core CPU performs enhanced correlation coefficient ECC mutual information estimation to obtain mutual information candidates; Merge phase-correlated candidates, geometric verification candidates, and mutual information candidates, and remove duplicates, retaining K candidates to form a candidate transform set {T_ij(1), T_ij(2), …, T_ij(K)}; The value of K ranges from 3 to 8.

3. The method according to claim 1, characterized in that, The method of constructing the potential energy ψ_ab(l_a, l_b) described in S4 further includes: For each edge (a, b) in the image pair graph E_pair, when a = (i, j) and b = (j, k) share the intermediate image j (i.e., forward sharing), the potential energy of the pair is ψ_ab(l_a, l_b) = λ · ‖T_ij(l_a) · T_jk(l_b) − T_ik*‖_F²; When a = (i, j) and b = (i, k) share the left diagram i, i.e., share in opposite directions, the potential energy is ψ_ab(l_a, l_b) = λ · ‖T_ij(l_a) · T_ik(l_b)⁻¹ − T_jk*‖_F²; The desired composite transformation T_ik* takes the candidate T_ik(argmin_l C_ik(l)) with the minimum pairwise matching cost in the candidate set of (i, k) when (i, k) ∈ E, i.e., i and k are adjacent image pairs; and takes the unit transformation I when (i, k) ∉ E. The value of λ ranges from 0.5 to 2.

0.

4. The method according to claim 1, characterized in that, The TRW-S discrete energy minimization solution described in S6 further includes: The MRF graph is decomposed into multiple subtrees, with each edge belonging to at least one subtree; Node traversal is divided into blocks based on image physical coordinates. Within each block, OpenMP multi-threaded parallelism is used, while between blocks, the topological order is followed to ensure monotonically decreasing energy. The K × K potential energy matrix for each edge is stored in row-major order, and the message update core loop is vectorized using Single Instruction Multiple Data (SIMD) instructions. The iteration stops when the relative energy change (E_(t-1) − E_t) / E_(t-1) between two adjacent iterations is less than the convergence threshold τ or the maximum number of iterations iter_max is reached, and the optimal label l_a* of each node is output.

5. The method according to claim 1, characterized in that, The abnormal node detection and candidate degradation regeneration described in S7 further include: For each node a, calculate the sum of potential energy under its optimal label S_a = Σ_(b ∈ N(a)) ψ_ab(l_a*,l_b*), where N(a) is the set of neighbors of node a in the image pair graph; When S_a is greater than the threshold θ_high, the node is marked as an abnormal node, where the value of θ_high ranges from 1.0 to 5.0; The candidate set of anomalous nodes is downgraded and regenerated, and the downgraded regeneration includes one or more of the following: expanding the SIFT matching radius, reducing the RANSAC threshold, and adding LOG polar coordinate candidates; Only the local iterations S3 to S6 are re-executed for the 2-hop neighborhood of the abnormal node to avoid recalculating the entire graph; If S_a is still greater than θ_high after local iteration, the image pair is deleted from the original image G and the corresponding position is set to zero in the sparse adjacency matrix Ã.

6. The method according to claim 1, characterized in that, The edge confidence vector w described in S8 is constructed as follows: w_a = exp(−S_a / σ²), where σ is a temperature parameter, ranging from 1.0 to 3.0; w_a has a dimension of 1 and a value range from 0 to 1. The closer w_a is to 0, the less reliable the candidate of the image pair is. The closer w_a is to 1, the more reliable the candidate of the image pair is. The edge confidence vector w serves as the initial robust weight for GNC annealing in the subsequent continuous domain global alignment module.

7. The method according to claim 1, characterized in that, The cascading execution of subsequent global alignment described in S9 further includes: S9.1 The sparse adjacency matrix Ã, the candidate transformation T_ij(l_a*) and the edge confidence vector w are delivered to the subsequent continuous domain global alignment module through the zero-copy shared memory of the Linux kernel DMA-BUF framework; S9.2 The global alignment module of the subsequent continuous domain executes the spectrum SE(2) synchronously with à and T_ij(l_a*) as input to obtain the initial value of the global pose; S9.3 The subsequent continuous domain global alignment module performs asymptotic nonconvex GNC annealing with w as the initial weight to further eliminate residual outliers in the continuous domain; S9.4 The final global pose is output by the subsequent continuous domain global alignment module and enters multi-band fusion.

8. The method according to claim 1, characterized in that, The heterogeneous resource scheduling described in S10 further includes: The small-core CPU performs tile image decoding, physical coordinate management, and global task scheduling; The GPU performs parallel computation of the FFT phase correlation and candidate generation stages using OpenCL. The large-core CPU performs MRF construction, TRW-S message passing, and abnormal node detection using OpenMP multi-threading. Tensor, image, and matrix data transfer between processing units is performed in a zero-copy manner through the Linux kernel DMA-BUF framework; The Dynamic Voltage and Frequency Scaling (DVFS) mechanism dynamically adjusts the voltage and frequency of each processing unit according to the current stage load: during the candidate generation stage, the GPU is fully loaded while the CPU is downclocked; during the TRW-S stage, the CPU is fully loaded while the GPU is downclocked.

9. The method according to claim 1, characterized in that, The method is implemented as a pre-discrete screening stage and a subsequent continuous domain global alignment module, wherein the subsequent continuous domain global alignment module includes a spectrum SE(2) synchronization module and a GNC annealing module; The sparse adjacency matrix Ã, candidate transformation T_ij(l_a*), and edge confidence vector w output by the method significantly expand the convergence domain of the subsequent continuous domain global alignment module and accelerate its convergence.

10. The method according to claim 1, characterized in that, The method is applied to three modalities: bright-field HE staining, immunohistochemical IHC staining, and fluorescence multichannel and Z-Stack three-dimensional tomography. The only difference between the different modalities is the replacement of the local feature sources generated as candidate samples in S1. Bright-field HE and IHC staining modalities were used as local feature sources by a combination of SIFT and normalized cross-correlation NCC. Fluorescence multichannel modalities are represented by a combination of SIFT and phase correlation as local feature sources; Z-Stack 3D tomographic mode uses a combination of SIFT and Sparse Microscopic Localization SML focal metric as a local feature source; The three modalities share the discrete candidate screening backends S2 to S9.

11. A discrete candidate screening device for pathological WSI images based on cycle consistency, characterized in that, include: The candidate generation module is used to acquire N tiles and generate K candidate relative transformations for each image pair; The image-pair graph construction module is used to construct Markov random fields (MRFs) using image pairs as nodes. The potential energy construction module is used to construct unary potential energy based on pairwise matching cost and pair potential energy based on cycle consistency residual; The discrete energy minimization module is used to solve for the optimal candidate label for the total energy of the MRF using the TRW-S algorithm; The abnormal node detection and degradation module is used to detect abnormal nodes and perform candidate degradation regeneration and local iteration on the node and its 2-hop neighborhood; The candidate selection result output module is used to output the sparse adjacency matrix Ã, the candidate transformation T_ij(l_a*), and the edge confidence vector w; The cascading execution module is used to deliver the candidate screening results to the subsequent contiguous domain global alignment module via shared memory; The heterogeneous scheduling module is used to allocate the tasks of the various modules to the neural network processing unit (NPU), graphics processing unit (GPU), large-core CPU, and small-core CPU for execution.

12. An electronic device comprising a processor and a memory, wherein the memory stores a computer program, characterized in that, When the processor executes the program, it implements the steps of the method according to any one of claims 1 to 10.

13. A computer-readable storage medium having a computer program stored thereon, characterized in that, When executed by a processor, the computer program implements the steps of the method according to any one of claims 1 to 10.

Images