A three-dimensional segmentation method for umbilical cord stem cells

Abstract

This invention relates to the field of biomedical image processing, and more particularly to a three-dimensional segmentation method for umbilical cord stem cells. The method includes: firstly, generating a continuous fluorescence density field based on the inversion of the imaging operator, and constructing a spatial metric tensor field using the point spread function and depth attenuation; secondly, determining the cell center and target mass under geodesic distance, constructing a power graph and optimizing the weights with semi-discrete optimal transfer to obtain mass-conserving three-dimensional units; subsequently, initializing a signed distance function with the three-dimensional units, and obtaining cell domains by combining data alignment, non-overlap, and topological constraints; finally, generating a synthetic density with the cell domains, using the imaging residual closed loop to drive division merging and metric updates, iterating until convergence, and outputting voxel-level masks and statistical indicators such as number, volume, spatial density, and adherent contact area.

Description

Technical Field

[0001] This invention relates to the field of biomedical image processing, and more particularly to a three-dimensional segmentation method for umbilical cord stem cells. Background Technology

[0002] Three-dimensional cell culture and mass production are becoming a key foundation for cell therapy and regenerative medicine. Production and quality control require voxel-level precision to obtain indicators such as cell number, volume, spatial density, adhesion contact area, and spatial distribution for batch-to-batch comparisons, expansion strategy optimization, and quality release. However, existing methods largely rely on two-dimensional slice inference or end-to-end segmentation networks, which struggle to handle anisotropic ambiguity caused by depth decay and point spread functions. They also struggle to maintain topological correctness in cell adhesion and dense contact scenarios and lack closed-loop verification consistent with imaging physics, leading to inconsistent statistical standards and insufficient batch-to-batch comparability. Summary of the Invention

[0003] To address the numerous problems existing in the prior art, this invention provides a three-dimensional segmentation method for umbilical cord stem cells. The method first obtains a continuous fluorescence density field through imaging operator inversion and constructs a spatial metric tensor field accordingly. Under geodesic distance, mass-conserving partitioning is achieved using a power graph and semi-discrete optimal transport. Then, a signed distance function is used in conjunction with data alignment, non-overlap, and topological constraints to obtain cell domains. Finally, imaging residual loop closure triggers cell division merging and metric updates. The method outputs voxel-level masks and standardized statistics to ensure physical consistency, topological correctness, and batch comparability.

[0004] A method for three-dimensional segmentation of umbilical cord stem cells includes the following steps:

[0005] A continuous fluorescence density field is generated by inverting the data of confocal imaging or light sheet imaging based on the imaging operator, and a spatial metric tensor field is constructed based on the point spread function and depth attenuation.

[0006] Under the geodesic distance induced by the spatial metric tensor field, the cell center and target mass are determined based on the continuous fluorescence density field. A power graph is constructed and the power graph weights are optimized through semi-discrete optimal transmission to obtain three-dimensional units and adjacency relationships that satisfy mass conservation.

[0007] The signed distance function is initialized with three-dimensional units under a spatial metric tensor field. Data alignment terms are constructed based on the continuous fluorescence density field and the signed distance function is jointly optimized. Non-overlapping constraints and topological constraints are applied to obtain cell domains and voxel-level masks.

[0008] Based on the synthetic density generated by the cell domain, synthetic observations are obtained through imaging operators and residuals are calculated. Based on the residuals, splitting, merging, or deformation is triggered and the power graph weights and spatial metric tensor field are updated. The process is iterated until convergence, and the final voxel-level mask is output.

[0009] Compared with the prior art, the advantages and beneficial effects of the present invention are as follows:

[0010] By inverting the imaging operator to generate a continuous fluorescence density field, explicit correction of depth attenuation and point spread function ambiguity is achieved, improving the separability of deep and axial details.

[0011] By constructing a spatial metric tensor field and geodesic distance using a point spread function and depth attenuation, an anisotropic geometric metric consistent with optical resolution is achieved, avoiding boundaries crossing microcarriers or low-confidence regions. By constructing a power graph under geodesic distance and performing semi-discrete optimal transmission weight optimization, a three-dimensional initial partition with single-cell mass conservation is achieved, overcoming volume distribution deviations caused by nearest-neighbor compression.

[0012] By jointly optimizing with a consistent signed distance function under the metric and imposing non-overlapping and topological constraints, mutually independent and simply connected cell domains are realized, reducing oversegmentation and undersegmentation;

[0013] By driving splitting, merging, and deformation based on residual closed-loop driven by imaging operators, and updating the power graph weights and spatial metric tensor field, physical consistency verification and adaptive correction are achieved, improving the interpretability and robustness of the results.

[0014] By unifying coordinates and measurements throughout the entire process and outputting voxels, standardized statistics of indicators such as quantity, volume, spatial density and wall contact area have been achieved, enhancing inter-batch comparability and the reliability of release decisions.

[0015] By using a two-channel design that separates and couples geometric constraints with mass conservation, robust segmentation of deep regions with weak signal-to-noise ratios is achieved, and dependence on specific network weights and data domain migrations is reduced. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the method of the present invention. Detailed Implementation

[0017] The embodiments of the present disclosure will now be described with reference to the accompanying drawings. However, it should be understood that these descriptions are exemplary only and are not intended to limit the scope of the disclosure. In the following detailed description, numerous specific details are set forth to provide a thorough understanding of the embodiments of the present disclosure for ease of explanation.

[0018] like Figure 1 As shown, a three-dimensional segmentation method for umbilical cord stem cells includes the following steps:

[0019] A continuous fluorescence density field is generated by inverting the data of confocal imaging or light sheet imaging based on the imaging operator, and a spatial metric tensor field is constructed based on the point spread function and depth attenuation.

[0020] Preferably, constructing the spatial metric tensor field includes forming an anisotropic metric based on the point spread function and depth attenuation, and introducing an uncertainty term generated by the residual and signal-to-noise ratio. When microcarriers are present, the motion cost is increased in the normal direction based on the distance field of the microcarrier surface to form an insurmountable constraint. Inverting to generate a continuous fluorescence density field includes performing flat field correction, stripe suppression and depth attenuation normalization on the volume data, and solving the non-negative regularized inverse problem based on the imaging operator.

[0021] Within the framework of the unification of imaging physics and geometry, this invention first uses volumetric data obtained from confocal imaging or light sheet imaging to generate a continuous fluorescence density field through inversion based on imaging operators. Simultaneously, it constructs a spatial metric tensor field based on the point spread function and depth attenuation. These two products serve as unified inputs for subsequent determination of cell centers, mass conservation partitioning, and surface optimization. This section explains the principles, implementation, and effects of this step and provides an example.

[0022] In the inversion stage, a correspondence is established between the original volume data and the optical response of the imaging process. A continuous fluorescence density field is obtained through a non-negative regularized variational solution, which recovers details while suppressing noise and fringe interference. The core solution objectives are as follows:

[0023]

[0024] Where ρ represents the continuous fluorescence density field, and y represents the volume data. The imaging operator is defined by θ, which represents the set of physical parameters such as the point spread function and depth attenuation. α and β are regularization weights, and φ is the penalty function for edge preservation. During the solution process, flat-field correction, fringe suppression, and depth attenuation normalization are performed first to reduce the influence of systematic inhomogeneities. Convergence is then achieved through layered multi-scale iterations. The technical advantage of this design lies in generating a density field with continuous intensity, clear boundaries, and consistent cross-layer distribution without altering the physical meaning, providing a reliable integrand for subsequent geodesic distance and mass volume integrals.

[0025] Spatial metric tensor fields are used to uniformly represent the anisotropy of point spread functions, resolution variations caused by depth decay, and uncertainties and geometric barriers derived from data quality in three-dimensional space. The synthesis of the metric follows the equation:

[0026]

[0027] Where G(x) is the spatial metric tensor field, x is the spatial position, and G... psf(x) is an anisotropy measure derived from the principal direction and principal axis length of the point spread function, where U(x) is the uncertainty field, I is the unit tensor, τ is the uncertainty weight, d(x) is the distance field at the microcarrier surface, n(x) is the microcarrier normal, and κ is the gain function that increases the motion cost along the normal. After the metric is defined, all distances, curvatures, and level set normal velocities are calculated under this metric. The technical effect is that it naturally penalizes paths traversing the microcarrier and low-confidence regions within the same geometry, maintaining anisotropy consistent with optical resolution.

[0028] The uncertainty field is used to reflect the reliability of data in the measurement. First, the generated continuous fluorescence density field is used to forward generate synthetic observations through imaging operators. Then, the difference between this and the volume data is calculated to obtain the residual. The residual is then fused with the signal-to-noise ratio estimated per channel to obtain the uncertainty field. This uncertainty only participates in the measurement modulation and does not interfere with the density field itself, thus avoiding bias in the inverted intensity.

[0029] The distance field on the microcarrier surface is obtained by threshold segmentation and surface reconstruction of the reflection or bright field channel, followed by calculation of the distance value and extraction of the normal using fast distance transformation. The smaller the distance, the larger the normal component, forming an impenetrable barrier. This avoids the segmentation boundary crossing the solid interface in cell adhesion or contact scenarios, ensuring the physical rationality of the spatial relationship.

[0030] The inversion process employs nonnegativity constraints and edge-preserving regularization to avoid negative intensity and over-smoothing. Gradient regularization suppresses noise, while second-order regularization maintains structural continuity. Since the imaging operator includes point spread functions and depth attenuation, the inversion results retain usable boundary contrast even at deeper layers. A multi-scale strategy improves convergence efficiency and avoids getting trapped in local minima.

[0031] The metric construction relies on two types of parameters: point spread function and depth attenuation. The point spread function parameters can be obtained by fitting local samples from fluorescent bead markers, or by using alternating minimization estimation in the absence of markers; the depth attenuation can be obtained by fitting the interlayer background trend. The resulting set of physical parameters remains unchanged in subsequent steps, with only the uncertainty weights updated using the residuals, without altering the fundamental physical properties of the point spread function and depth attenuation.

[0032] The comprehensive effects of the above design are reflected in three aspects. First, the continuous fluorescence density field provides a stable intensity basis for determining the cell center and target mass, with a one-to-one correspondence between the mass volume component and the density field. Second, the spatial metric tensor field supports geodesic distance, power graph generation, and level set advancement with unified geometric constraints, significantly reducing the probability of missegmentation in wall-hugging scenes and low-confidence regions. Third, residual-driven uncertainty modulation forms a closed-loop mapping from data quality to geometric cost, making subsequent steps more robust to complex structures.

[0033] Example: A batch of confocal volumetric data was processed. First, flat-field correction and fringe suppression were performed, followed by attenuation normalization along the depth direction to obtain preprocessed volumetric data. An imaging operator was constructed using the preprocessed volumetric data, the calibrated point spread function, and depth attenuation parameters. Inversion was performed according to the aforementioned objective function to obtain a continuous fluorescence density field. The continuous fluorescence density field was forward-generated into a synthetic observation using the imaging operator, and the difference between this observation and the volumetric data was used to obtain the residual. This residual was then combined with the signal-to-noise ratio obtained through layer-by-layer statistics to form an uncertainty field. Using imaging data of fluorescent microspheres at different depths and resolutions, a three-dimensional model was used to estimate the principal axis direction and length of the point spread function, and an anisotropy metric component was formed accordingly. Threshold segmentation was performed on the bright-field channels of the same batch, and the microcarrier surface was reconstructed. The distance field and normal were calculated. The anisotropy metric, uncertainty component, and microcarrier normal barrier were combined according to the aforementioned synthetic formula to obtain a spatial metric tensor field. At this point, the continuous fluorescence density field and spatial metric tensor field have been generated, serving as unified inputs for subsequent determination of cell center and target mass, construction of a power graph under geodesic distance, mass conservation partitioning through semi-discrete optimal transport, and optimization of a consistent signed distance function under metric. In this embodiment, inversion and metric construction are performed within the same coordinate system to avoid cross-scale interpolation errors; in the adherent region, the metric normal cost is significantly increased, and subsequent surface evolution will not traverse the microcarrier, thus maintaining the rationality of topological and contact relationships.

[0034] Under the geodesic distance induced by the spatial metric tensor field, the cell center and target mass are determined based on the continuous fluorescence density field. A power graph is constructed and the power graph weights are optimized through semi-discrete optimal transmission to obtain three-dimensional units and adjacency relationships that satisfy mass conservation.

[0035] This invention transforms a continuous fluorescence density field into a three-dimensional cell partition that satisfies mass conservation under geodesic distances induced by a spatial metric tensor field. The core process involves cell center determination, target mass setting, geodesic distance calculation, power graph construction, and semi-discrete optimal transport weight optimization, outputting the three-dimensional cells and their adjacency relationships. This process maintains the same coordinates and geometric metrics as the previously generated continuous fluorescence density field and spatial metric tensor field, avoiding cross-domain inconsistencies.

[0036] Cell centers were determined based on nuclear channels within a continuous fluorescence density field. A superlevel set sequence was constructed, and topological persistence was calculated. Local maxima with high persistence were selected as candidate points. To eliminate overly dense candidate points, non-maximum suppression and nearest-neighbor merging under geodesic distances induced by the spatial metric tensor field were employed to obtain the set of cell centers. This approach maintains a reasonable minimum spacing at points of metric principal axis stretching or compression, thus ensuring geometric consistency with subsequent geodesic distance-based partitioning.

[0037] The target mass is used to constrain the mass conservation of the three-dimensional unit. Taking each cell center as the core, the volume integral of the cytoplasmic channels within the continuous fluorescence density field in the Earth's neighborhood is performed to obtain an initial mass value. This initial mass value is then constrained by combining the value range obtained from batch statistics, forming a target mass set. If cytoplasmic channels are lacking, the intensity integral of the available channels is linearly mapped to the initial mass value based on the calibrated channel mapping before further constraining. Through this strategy, the source of the target mass and subsequent mass constraints maintain traceable consistency.

[0038] The geodesic distance is calculated within a spatial metric tensor field, using an anisotropic fast travel method to obtain the geodesic distance field from each cell center to the global domain. Boundaries are set with unreachable or high-cost conditions to avoid irrelevant out-of-bounds travel. This distance propagates along low-cost directions in regions dominated by point spread function anisotropy and depth attenuation, thus providing geometric guidance consistent with the imaging resolution for the subsequent boundary positions.

[0039] The power graph is constructed based on the geodesic distance field and power graph weights, and voxel-level partitioning is performed using the power distance function:

[0040] ∏ i (x)=D i (x) 2 -w i

[0041] Among them, ∏ i (x) is the power-law distance function of position x with respect to the center of the i-th cell, D i (x) represents the geodesic distance under the spatial metric tensor field, w i The weights are determined by the power graph. For each location, the index of the minimum power distance is used to form a set of 3D cells, and adjacency relationships are extracted on the voxel label boundaries. This construction method is geometrically equivalent to introducing an adjustable bias into a potential field of squared geodesic distances, preserving spatial anisotropy while providing degrees of freedom for subsequent mass matching.

[0042] Semi-discrete optimal transport is used to optimize power graph weights to match target mass and minimize geodesic displacement cost under fixed cell centers and spatial metric tensor fields. The optimization objective and constraints are:

[0043]

[0044] Among them, L i (w) represents the weights of the power graph {w} i The i-th three-dimensional unit is determined, where ρ(x) is a cytoplasmic channel of continuous fluorescence density field, and m i For the target quality, the gradient of the power graph weights is:

[0045]

[0046] in, To optimize the objective, Let be the mass within the current 3D element. Iteration is performed using Newton's method or a quasi-Newton method, with the 3D element and mass recalculated after each update, until the mass deviation meets the preset tolerance and the objective function no longer decreases significantly. In the above expression, x is the spatial position, ρ(x) is the intensity of the continuous fluorescence density field, and D... i (x) represents the geodesic distance, w i For the power graph weights, m i For target quality, M i (w) represents the current quality. To optimize the objective.

[0047] After obtaining the 3D elements, the adjacency graph can be directly constructed based on the common boundaries of the voxel labels, and the positions of triangular and multi-pronged contacts can be detected. This adjacency relationship is used in subsequent surface optimization to construct non-overlapping constraints and multi-body repulsion between adjacent elements, ensuring local consistency in boundary advancement.

[0048] The technical benefits of the above process are reflected in three aspects. First, the cell center obtains a stable position under the dual constraints of topology and geodesic distance, reducing unnecessary subsequent splitting and merging. Second, the power graph is generated under the spatial metric tensor field, and the boundary automatically fits the geometric environment after resolution anisotropy and uncertainty modulation. Third, the semi-discrete optimal transmission uses mass conservation as a hard constraint, while minimizing the geodesic displacement cost, making the three-dimensional unit correspond one-to-one with the target mass, providing a consistent input for the initialization of the subsequent signed distance function and mass bias.

[0049] Example: For a set of confocal data, the continuous fluorescence density field and spatial metric tensor field generated in the previous steps are first used. A hyperlevel set sequence is constructed on the nuclear channels, and topological persistent cohomology is calculated. Local maxima with high persistence are selected as candidates. Subsequently, non-maximum suppression and nearest neighbor merging are performed under geodesic distance to obtain a set of cell centers. Under geodesic distance induced by the spatial metric tensor field, the set of cell centers is determined. i Centered on, radius r i Earth B measurement G (s i ,r i )={x∣d G (x,s i )≤r iThe cytoplasmic channels of the continuous fluorescence density field are volume-integrated, and a target mass set is formed with reference to historical batch intervals. Anisotropic fast traversal is used to calculate the geodesic distance field from the center of each cell to the entire region. The power distance function is constructed by subtracting the power graph weights from the squared geodesic distance. Three-dimensional cells are divided according to the minimum power distance, and adjacency relationships are extracted. With mass conservation as a constraint, the geodesic displacement cost is minimized, and the power graph weights are iteratively optimized. After convergence, a three-dimensional cell set and adjacency graph consistent with the target mass are obtained. In this embodiment, the cell boundaries in the adherent region bend along a low-cost direction without crossing geometric barriers, maintaining a reasonable cell volume distribution in deeper locations with weaker signals, thus meeting the input requirements for subsequent surface initialization and joint optimization based on the signed distance function.

[0050] Preferably, determining the cell center involves constructing a superlevel set sequence on the nuclear channel of a continuous fluorescence density field, screening local maxima as candidate points based on topological persistent homology, and performing nonmaximum suppression and nearest neighbor merging under geodesic distance induced by a spatial metric tensor field to obtain the cell center.

[0051] This invention determines the set of cell centers based on the nuclear channels of a continuous fluorescence density field under geodesic distance induced by a spatial metric tensor field, and outputs this set for subsequent quality setting, power graph construction, and surface initialization. The core of this step lies in determining the stable maxima locations using topological persistent homology and performing nonmaxima suppression and nearest neighbor merging using geodesic distance, ensuring that the cell centers maintain geometric consistency with the subsequent geodesic distance-based region partitioning.

[0052] First, scale-consistent peak enhancement is performed on the nuclear channels of the continuous fluorescence density field as a preprocessing step for candidate detection, without altering the definitions of subsequent metrics and distances. Then, a sequence of intensity hyperlevel sets is constructed, and zero-dimensional topological persistent cohomology is used to track the generation and merging of connected components. Each local maximum is associated with a pair of birth and death thresholds, and the difference between the two is used as persistence. To avoid lengthy explanations of well-known theories, this invention uses persistence only at the application level as a quantitative indicator of maximum stability to eliminate spurious peaks caused by noise or fringes.

[0053] To facilitate reproduction, the core computational expressions used in this invention are presented. The intensity of the nuclear channel is denoted as I(x), where x is the spatial location. The hyperlevel set K is defined. τ ={x|I(x)≥τ}, for the corresponding local maxima in topological filtering, let the birth threshold be τ. b The extinction threshold is τ d Then persistence is defined as p = τ b -τ dWhere p is the persistence. The candidate score q = λp + (1-λ)I(x) is composed of persistence and peak intensity, where q is the candidate score and λ is the weighting coefficient. The spatial metric tensor field is denoted as G(x), and the geodesic distance is defined as:

[0054]

[0055] Where γ is a differentiable curve connecting x and y. This represents the curve velocity.

[0056] After obtaining candidate points and their corresponding candidate scores, metric-consistent nonmaximum suppression is performed. The criterion is: if the geodesic distance between two candidate points is less than the set minimum geodesic interval, the candidate with the higher score is retained, and the other is deleted. This operation is also effective in regions with significant anisotropic resolution and depth attenuation, because the geodesic distance has already absorbed the anisotropy and uncertainty cost of the spatial metric tensor field. Next, nearest neighbor merging is performed: hierarchical clustering is constructed using geodesic distance as the metric, the cluster radius is a multiple of the minimum geodesic interval, and the representative position of the cluster is the geometric centroid weighted by the candidate scores, thereby integrating local multimodalities into a single cell center and avoiding abnormally fragmented units in the subsequent power graph.

[0057] The above process maintains consistency with the preceding steps in terms of coordinates and metrics, avoiding deviations caused by cross-scale or cross-coordinate resampling. Topological persistent cohomology provides stable extreme value screening for threshold changes, while geodesic distance-dominated suppression and merging ensures intrinsic consistency between the center interval and subsequent geodesic distance-based region segmentation. The combination of these two approaches gives the cell center set both topological stability and geometric consistency, reducing false positives and false negatives under conditions of weak signals, adherence to walls, or significant anisotropy.

[0058] In exceptional circumstances, if a candidate point with high persistence appears but is located in a strong stripe residue or saturation region, it can be suppressed by introducing a confidence reduction weight from the prior residual and signal-to-noise ratio into the candidate score. This reduction weight does not change the definition of persistence; it only lowers the priority of the candidate point during the sorting stage, maintaining the physical consistency and interpretability of the algorithm.

[0059] The output of this invention is a set of cell centers and their candidate scores. The data structure includes two parts: position and score. The position is given by three-dimensional coordinates, and the score is used for subsequent adaptive selection of the neighborhood scale for target quality and backtracking of anomalous units. The output results are directly used in the next step of solving the geodesic distance field and constructing the power graph, avoiding redundant calculations and information loss.

[0060] Example: A set of samples containing both adherent and suspended cells was selected. Using the previously obtained continuous fluorescence density field and spatial metric tensor field, scale-consistent peak enhancement was first performed on the cell nuclear channels. Subsequently, an intensity hyperlevel set sequence was constructed, and zero-dimensional topological persistence cohomology was calculated to obtain the persistence of each local maximum. The persistence and peak intensity were linearly fused to form a candidate score, and all candidate points were sorted in descending order of candidate score. Geodesic distances were calculated in the spatial metric tensor field using the anisotropic fast travel method, and non-maximum suppression was performed on the sorted list, with the minimum geodesic interval set according to the physical scale. On the remaining candidate point set, hierarchical merging was performed based on the geodesic distance, with the cluster radius on the same order of magnitude as the minimum geodesic interval, and the representative point position was taken as the geometric centroid weighted by the candidate score. The final set of cell centers covered the high-confidence positions of the cell nuclei, maintained a reasonable relationship between the center spacing and the microcarrier geometry in the adherent region, and retained stable centers in the deep low signal-to-noise ratio region. By directly passing this set to the subsequent target quality setting and power graph construction, the number of subsequent splitting and merging operations can be reduced, and the geometric assumptions and physical constraints consistent with the spatial metric tensor field can be maintained throughout the overall process.

[0061] Preferably, determining the target quality includes performing volume integral of the cytoplasmic channels of the continuous fluorescence density field within the geodesic neighborhood induced by the spatial metric tensor field with the cell center as the core, and limiting the amplitude according to the value range determined by historical batch statistics. In the absence of cytoplasmic channels, the target quality is determined by the intensity integral of the available channels based on the calibrated channel mapping.

[0062] This invention uses geodesic distance induced by a spatial metric tensor field as a unified geometry. Centering on the cell center, it performs volume integration on the cytoplasmic channels of a continuous fluorescence density field to obtain the target mass for subsequent mass conservation partitioning. The overall approach is as follows: In a physical and geometric coordinate system consistent with the preceding steps, a geodesic neighborhood is defined. The voxel intensities related to the cytoplasm within this neighborhood are extracted as initial mass values, and then amplitude is limited by combining statistical intervals from historical batches. When cytoplasmic channels are lacking, the target mass is obtained by integrating the intensities of available channels through calibrated channel mapping. This design ensures that the target mass possesses both local optical consistency and comparability with global statistics across process batches.

[0063] The core calculation uses the following expression. Let s be the cell center. i Let the spatial metric tensor field be G(x); and let the geodesic distance under this metric be d. G (x,s i ); using cytoplasmic channels of continuous fluorescence density field as ρ cyt (x); with the voxel physical volume as ΔV; and the geodesic neighborhood radius as r. i The initial mass value is:

[0064]

[0065] Where x is the spatial position, s i Let G(x) be the spatial metric tensor field at the i-th cell center, and d be the metric tensor field. G (x,s i ) represents the geodetic distance, ρ cyt (x) represents the cytoplasmic channel strength, ΔV represents the voxel physical volume, and r i To measure the radius of the geodetic neighborhood, This is the initial mass value. The geodesic neighborhood radius is adaptively correlated with the local scale, using r... i =γr nuc,i Where γ is a dimensionless proportionality constant, r nuc,i Local scale estimation of the cell nucleus is performed. The local scale can be obtained from the distance transformation of the nuclear channels or the Laplacian-Gaussian response; its general implementation is not detailed here. To mitigate batch-to-batch imaging gain, staining differences, and background shift, linear correction and statistical limiting are used to obtain the target quality.

[0066]

[0067] Where a and b are linear correction coefficients, and [L,U] is the value range determined by historical batch statistics, clip(z,[L,U])=min(max(z,L),U). The linear correction is obtained by global regression of the same batch of data, and the value range is determined by the quantile statistics of historical qualified batches, which avoids the interference of outliers on subsequent constraints and maintains comparability across batches.

[0068] In the absence of cytoplasmic channels, the target mass is obtained through calibrated channel mapping. Let the available channel strength be ρ. alt (x); with mapping coefficients α and β; then we have:

[0069]

[0070] Where, ρ alt (x) can represent the intensity of nuclear channels or cell membrane channels. α and β are obtained by fitting historical synchronous calibration data, maintaining consistent scales between channels. Subsequently, linear correction and statistical limiting are performed according to the previous formula to obtain m. i When the available channels are cell membrane channels, shell weights can be assigned to locations near the cell membrane within the geodesic neighborhood to improve sensitivity to boundary signals; this weight only affects The formation of this does not change the limiting and correction process.

[0071] The geodesic neighborhood defined using a spatial metric tensor field has two direct effects. First, the neighborhood shape automatically adapts to the anisotropy of the point spread function and variations in depth resolution, avoiding the introduction of excessive irrelevant voxels when axial resolution is low. Second, in the presence of microcarriers or low-confidence regions, the geodesic distance increases the passage cost, causing the neighborhood to naturally shrink to a reliable region, reducing the quality bias caused by spurious signals. This is geometrically consistent with subsequent power graph and surface optimization based on the same metric, avoiding inconsistencies in scale and shape across steps.

[0072] To improve robustness, one can... A multi-scale robust estimation method is employed. Let the radius sequence be... The corresponding initial value set is Then take:

[0073]

[0074] in, We take a finite number of nearest neighbor values ​​in the local scale, and median(·) is the median operator. This process is more stable when there are fluctuations in boundary intensity or local occlusion. The above robust steps do not change the overall process, only the way the initial values ​​are calculated is replaced.

[0075] The obtained target mass directly serves as a hard constraint for subsequent semi-discrete optimal transport, providing a definite mass quota for the optimization of the power graph weights. Since the target mass originates from the geodesic neighborhood volume integral consistent with subsequent metrics, the adjustment of weights during power graph iteration will not deviate from a physically reasonable mass scale, avoiding overexpansion or collapse.

[0076] Example: On a set of confocal data containing microcarriers, a continuous fluorescence density field and a spatial metric tensor field are first obtained from the preprocessor sequence, followed by a cell center step to obtain a set of cell centers. Using each cell center as a core, the geodesic distance field is calculated using the anisotropic fast travel method, and a geodesic neighborhood is constructed according to an adaptive radius. Within the geodesic neighborhood, volume integrals are performed on the cytoplasmic channels to obtain initial mass values. Subsequently, a value interval is set based on the quantile statistics of historical qualified batches, and linear correction and amplitude limiting are performed to obtain the target mass. When individual samples lack cytoplasmic channels, nuclear channels are used as available channels. Mapping coefficients are obtained through intra-batch regression, and initial mass values ​​are calculated and amplitude limited. The final target mass set and the cell center set are then fed into the power graph construction and semi-discrete optimal transmission steps. The power graph weights achieve mass conservation after a few iterations, and the adjacency relationships are stable. Subsequent surface initialization and joint optimization processes no longer require significant mass compensation. The above embodiments show that the combination of geodesic neighborhood volume integration and statistical limiting significantly reduces the quality deviation caused by imaging gain, local occlusion and attached structures, ensuring the reproducibility of the segmentation process in different batches and under different culture conditions.

[0077] Preferably, constructing a power graph involves solving for the geodesic distance field to the center of each cell under a spatial metric tensor field, constructing a power distance function with the square of the geodesic distance and the power graph weights, generating a power graph by dividing it according to the minimum power distance, and obtaining the adjacency relationship.

[0078] This invention constructs a power graph under geodesic distance induced by a spatial metric tensor field, completing the voxel-level 3D cell partitioning and outputting adjacency relationships. The goal is to control the voxel assignment of the control unit in physical and geometric coordinates consistent with the previous process using adjustable power graph weights, providing stable and traceable initial values ​​for subsequent mass conservation and surface initialization.

[0079] Given the set of cell centers, the global geodesic distance field is first solved for each cell center under a spatial metric tensor field. To ensure consistency with the anisotropic resolution and depth attenuation of the previous step, anisotropic fast travel is employed. During the calculation, the microcarrier normal barrier and uncertainty cost are used as intrinsic metrics, without introducing additional geometric assumptions inconsistent with the previous step. Discrete implementation preserves the original voxel spacing and avoids isovoxel resampling to prevent quadratic interpolation errors in distance and mass. Boundary conditions are set as unreachable or high-cost regions outside the bounding box to prevent curves from crossing the sampling domain.

[0080] The power graph uses the power distance function, which is the square of the geodesic distance and the power graph weights, as the criterion for determination, and the voxel assignment adopts the minimum power distance principle. The core expression is as follows:

[0081] ∏ i (x)=D i (x) 2 -w i

[0082]

[0083] Among them, ∏ i (x) is the power-law distance function of position x with respect to the center of the i-th cell, D i (x) represents the geodesic distance from the center of the i-th cell to position x under the spatial metric tensor field, w i Let L be the weight of the i-th power graph. i (w) represents the i-th 3D unit determined by the weight vector w. The voxel labels are... It is determined that when there is a tie for the minimum value, the one with the smaller geodesic distance is given priority. If they are still tied, the index is used to break the tie and ensure consistency in repeated runs.

[0084] The initialization of the power graph weights follows the principle of compatibility with subsequent mass conservation. If the target mass has been obtained through previous steps, the weights are given initial values ​​proportional to the logarithm of the target mass; if only geometric initialization is performed, the weights are zero vectors. The initial values ​​do not change the form of the power graph, but only affect the subtle assignment of voxels among several cell centers, which will be updated iteratively through semi-discrete optimal transport.

[0085] To obtain stable adjacency relationships, the voxel label field is scanned once. For adjacent voxels with label changes, the adjacent cell center index pairs are recorded, and duplicates are deduplicated using a hash method to form an adjacency list. Simultaneously, at 3D intersections, triangular and multi-branch contacts are statistically analyzed and their spatial positions are recorded. For subsequent surface initialization and normal estimation, on the common boundary voxel set of each adjacent index pair, the empirical mean of the local normal and a rough estimate of the curvature are calculated. These geometric statistics do not change the labels and are only passed as priors to subsequent steps.

[0086] To suppress isolated small fragments caused by numerical noise, consistency cleaning is performed after constructing the power graph. Specifically, connected component labeling is performed on each 3D cell, and fragments with areas below a threshold are merged into their nearest neighboring cells. For extremely thin and elongated fragments, a voxel-level morphological correction is performed along the normal direction of the common boundary, assigning them to neighboring cells with lower normal energy. This cleaning does not change the power graph's decision criteria; it only stabilizes areas of stagnation and numerical fluctuation.

[0087] The above design yields three direct effects. First, geodesic distance is defined under the spatial metric tensor field, and the cell boundaries naturally conform to anisotropic resolution and depth attenuation, avoiding unreasonable cross-layer connections when axial resolution is insufficient. Second, the power graph weights provide adjustable biases, enabling mass conservation through weight iteration while maintaining geometric compatibility, avoiding systematic volume deviations caused by relying on simple nearest distance. Third, voxel-level adjacency relationships and the common boundary voxel set provide structural indexes for subsequent non-overlapping constraints and multi-body repulsion, reducing topological oscillations during the surface optimization stage.

[0088] Example: On a set of volume data containing a mixture of adherent and suspended cells, the geodesic distance field of all cell centers is first calculated using the anisotropic fast travel method. Using the zero vector as the initial value of the power graph weights, voxels are partitioned according to the power distance function to generate initial 3D units. Subsequently, the voxel label boundaries are scanned, an adjacency list is constructed, and the triangular contact positions are recorded. Simultaneously, the set of voxels representing the common boundary is saved for subsequent coarse estimation of normals and curvature. Connectivity labeling is performed on each 3D unit, removing small fragments and incorporating lower-energy neighboring units. In this example, the boundary of the adherent region bends along the high-cost direction of the microcarrier normal rather than crossing it; the units in the suspended region are approximately spherical and continuous between layers; in locations with low signal-to-noise ratios at depth, the power graph boundary coincides with the low-cost valley of the geodesic distance, and no long, thin channels crossing low-confidence regions appear. The final 3D units and adjacency relationships are directly passed to the semi-discrete optimal transport step to update the power graph weights and used as geometric input for surface initialization.

[0089] Preferably, the optimization of semi-discrete optimal transport includes performing convex optimization iteration on the power graph weights, so that the volume fraction of the continuous fluorescence density field in the cytoplasmic channel within each three-dimensional cell is equal to the corresponding target mass, and generating a mass-conserving dual variable for the mass bias of the signed distance function.

[0090] This invention, under the joint constraints of the previously generated continuous fluorescence density field and the spatial metric tensor field, performs convex optimization iterations on the power graph weights to ensure that the cell mass within each 3D cell is consistent with the target mass, and outputs a mass-conserving dual variable for subsequent signed distance functions to apply mass bias. This step uses the geodesic distance induced by the spatial metric tensor field as the geometric basis, maintaining consistent coordinates and metrics with the power graph construction and surface optimization, without changing the already determined cell centers and adjacency relationships.

[0091] Power graphs use a power distance function to determine voxel assignment, defined as follows:

[0092] ∏ i (x)=D i (x) 2 -w i

[0093] Among them, ∏ i (x) is the power-law distance function of position x with respect to the center of the i-th cell, D i (x) represents the geodesic distance from the center of the i-th cell to position x under the spatial metric tensor field, w i The weights are power-law graph weights. The set of three-dimensional cells L is obtained by partitioning using the minimum power-law distance. i (w). Given the power graph weights, the mass of the cytoplasmic channel in the i-th three-dimensional unit of the continuous fluorescence density field is:

[0094]

[0095] Where, ρ cyt (x) represents the cytoplasmic channel strength. The goal is to achieve mass conservation with minimal geodesic displacement cost. The optimization problem is denoted as:

[0096]

[0097] Where, m i The target quality obtained from the previous section; M i (w)=m i The semi-discrete optimal transport, under this formulation, can be transformed into a convex optimization of power graph weights, whose first derivative has a closed form:

[0098]

[0099] in, Let be the objective function. This formula shows that the direction of each weight update is directly given by the "deviation between the current quality and the target quality", ensuring the physical interpretability of the iteration process.

[0100] To improve convergence speed and control numerical stability, second-order information from the sparse approximation is introduced. The common boundary of adjacent 3D element pairs is discretized, and the non-zero coupling terms of the Hersey matrix are obtained using the standard approximation of boundary integrals. Second-order terms of non-adjacent element pairs are ignored. Based on this, Newton's method or a quasi-Newton method is used for iteration, combined with backtracking search to ensure that the objective function monotonically decreases. Numerical estimation of quality and gradient is achieved by voxel summation and incremental updates of the boundary voxel set, avoiding global recalculation in each iteration. For empty elements, the corresponding power graph weights are reduced and a small number of voxels are recovered by referring to adjacency relationships, and then iteration continues. For elements with consistently high quality and multiple stable peaks within the element, they are marked as needing splitting and handled in subsequent steps.

[0101] To link with the mass bias voltage of the subsequent signed distance function, a mass conservation dual variable is introduced and updated in an augmented form:

[0102] λ i ←λ i +η(M i (w)-m i )

[0103] Where, λ i Let η be the mass conservation dual variable for the i-th unit, and η be the dual step size. The dual variable does not change the convex optimization trajectory of the power graph weights, but rather acts as an external bias for the normal velocity during the surface evolution stage, prompting the geometric boundary and mass constraints to reach consistency at local details.

[0104] At the implementation level, voxel assignment is given by the minimum index of the power-law distance function; the mass integral is implemented through voxel parallel summation; the first-order gradient is directly given by the mass deviation of the 3D cells; the second-order approximation is calculated only on cell pairs indexed by the adjacency list, and its complexity is linearly related to the adjacency degree. The termination condition of the line search adopts a dual standard of constraint on the relative reduction of the objective function and the mass error to prevent label oscillations caused by large step sizes. In terms of numerical accuracy, all voxel operations retain the voxel physical volume weights to ensure that the anisotropy generated by the voxel mesh does not distort the mass statistics.

[0105] This step yields three benefits. First, mass conservation is precisely achieved at the power graph weight level, providing stable boundary conditions for subsequent geometric evolution and reducing the sensitivity of the geometric stage to intensity noise. Second, the minimum geodesic displacement cost allows for the redistribution of mass in the "shortest" manner under the spatial metric tensor field, avoiding crossing microcarriers or low-confidence regions. Third, the separation of dual variables and weight optimization allows "geometric adjustment" and "mass constraints" to form complementary pathways, improving the reliability and interpretability of overall convergence.

[0106] Example: On samples containing adherent and suspended cells, the cell centers, target mass, and initial values ​​of the power graph obtained in the previous step are used. The power graph weights are initialized with zero vectors, and the voxel assignments and 3D cell masses are calculated to form the first-order gradient. The sparse coupling term of the Hersey matrix is ​​estimated on the common boundary voxel set using an adjacency list, and Newton iteration is performed in conjunction with backtracking search. After each update, only the affected 3D cells and their adjacency pairs are incrementally redistributed and their mass is statistically analyzed. After approximately 10 to 30 iterations, the mass deviation of all 3D cells enters the preset tolerance, the relative decrease of the objective function is less than the threshold, and the dual variable is stable. No mass migration across the microcarrier occurs in the adherent region, and the volume distribution of 3D cells in deep, low signal-to-noise ratio locations is consistent with the target mass. The final output power graph weights, 3D cells, and mass-conserving dual variables are directly used for surface initialization and mass bias of the signed distance function. Subsequent geometric stages only require local fine-tuning to obtain cell boundaries consistent with imaging physics.

[0107] The signed distance function is initialized with three-dimensional units under a spatial metric tensor field. Data alignment terms are constructed based on the continuous fluorescence density field and the signed distance function is jointly optimized. Non-overlapping constraints and topological constraints are applied to obtain cell domains and voxel-level masks.

[0108] This invention uses three-dimensional units as geometric initial values, initializes a signed distance function under a spatial metric tensor field, constructs a data alignment term based on a continuous fluorescence density field, and jointly optimizes it with geometric regularization, non-overlapping constraints, and topological constraints to finally obtain the cellular domain and voxel-level mask. All calculations are performed in coordinates and metrics consistent with the previous steps, ensuring consistency in geodesic distance, quality statistics, and boundary evolution.

[0109] Let the set of voxels of the i-th 3D unit be the initial interior domain, and let the signed distance function be f. i (x) where x is the spatial location. Initialization is obtained by solving the anisotropic distance equation under the spatial metric tensor field: The inner domain takes negative values, the outer domain takes positive values, and a narrow band is established around the zero level set. Where G(x) is the spatial metric tensor field, and w is the narrowband thickness. The cellular domain Ω is thus defined. i ={x|f i (x)≤0} and normal direction and the mean curvature κ calculated consistently under the metric i (x).

[0110] Data alignment is implemented in two categories based on available channels. When cell membrane channels exist, the cell membrane probability map is denoted as P. mem (x), using boundary attraction energy:

[0111]

[0112] Where φ is a convex function. When only cytoplasmic channels exist, a region-consistent energy is used:

[0113]

[0114] Where, ρ cyt (x) represents a cytoplasmic channel with a continuous fluorescence density field. in With l out This is the log-likelihood term. The two types of data alignment terms are chosen individually or in a weighted combination in the implementation. Geometric regularization includes distance function regularization and curvature regularization:

[0115]

[0116] Non-overlapping constraints are implemented with soft penalties, denoted by the smooth penalty function σ(·):

[0117]

[0118] The term is evaluated on the index pairs given the adjacency relation to suppress voxel overlap between 3D units.

[0119] Mass conservation provides bias pressure on the boundary through the mass conservation dual variable. Let the mass of the i-th element under the current geometry be:

[0120]

[0121] The target mass is m i The mass bias velocity term is defined as v mass,i (x)=β(M i (f)-m iWhere β is the weight; the mass-conserving dual variable is denoted as λ. i It also performs rolling updates to transmit quality deviation information between multiple iterations.

[0122] To avoid tangential perturbations and obtain numerically stable normal propagation, a level set evolution form is adopted:

[0123]

[0124] Where, δ τ For a smoothed Dirac function, the synthesized normal velocity is:

[0125]

[0126] v data,i α is obtained from the shape derivative of the data alignment term. dist With α curv As the weight. Every few steps, f is adjusted. i Perform a reinitialization based on the fast-moving method, so that Keep it close to 1 and update the narrowband and geometry in sync.

[0127] Topological constraints aim for simple connectivity and absence of holes. Within a narrow band, changes in the Eulerian characteristic are monitored, and unwanted splitting, merging, or puncturing events are prohibited. Once a risky event is detected, the system reverts to the most recently consistent state and increases the weights of curvature regularization and non-overlap penalties until the event disappears. This process does not change the energy form, only adjusts the weights and step size, ensuring topological stability of the evolution.

[0128] Once all signed distance functions converge, output the set of cell domains Ω. i With voxel-level mask Boundary uncertainties and event logs are retained for use in subsequent steps. The direct effects of this step are: the boundary is consistent with the cell membrane signal or cytoplasmic intensity model; there is no overlap between cell domains and they remain simply connected; mass bias promotes geometric consistency with the target mass; and under the guidance of the spatial metric tensor field, the boundary does not cross microcarriers or low-confidence regions.

[0129] Example: Using volumetric data containing both adherent and suspended regions, the previously obtained 3D cells, adjacency relationships, target mass, and mass-conserving dual variables are employed. First, the distance equation for each 3D cell is solved under a spatial metric tensor field to obtain an initial signed distance function and establish a narrow band with a thickness of several voxels. For data alignment, boundary attraction energy is used when cell membrane channels exist, and region consistency energy is used when only cytoplasmic channels exist; the weights for both types are set based on channel mass. Subsequently, the level set is advanced with a synthetic normal velocity, re-initialized every 10 to 20 steps, and non-overlap penalties are evaluated on adjacency index pairs; potential topological events are resolved through backtracking and weighting. After approximately several tens of iterations, the cell domain boundary stabilizes near the high-probability band of the cell membrane, the boundary of the adherent region bends along the high-cost direction of the microcarrier normal but does not cross it, and smooth, volume-integral-consistent boundaries are obtained in deep, weak signal areas using mass bias and curvature regularization. Finally, a voxel-level mask and cell domain set are derived, providing geometric input for subsequent imaging consistency closure and statistical index calculations.

[0130] Preferably, initializing the signed distance function involves solving the Eikonal equation for distance calculation at the boundary of the three-dimensional element to obtain the initial signed distance function, and updating the level set within a narrow band of fixed thickness.

[0131] This invention uses three-dimensional units as the initial geometric domain, constructs a signed distance function under a spatial metric tensor field, and updates the level set within a narrow band of fixed thickness, thereby providing a numerically stable and physically consistent boundary representation for subsequent joint optimization and constraint application. This step uses the same coordinates and metric as the previously generated continuous fluorescence density field and spatial metric tensor field, avoiding cross-domain errors.

[0132] Let the set of voxels of the i-th 3D unit be denoted as the initial interior region, and its boundary be the initial zero level set. The signed distance function is denoted as f. i (x), where x is the spatial location. Under the metric of the spatial metric tensor field denoted as G(x), the initial f i (x) is defined by the anisotropic Eikonal equation: At the initial boundary, in the formula, For spatial gradient, G -1 (x) represents the inverse of the spatial metric tensor field. The inner domain is negative and the outer domain is positive. Numerically, an anisotropic fast-marching method is used for solution, and gradient discretization uses unidirectional upwind differencing to ensure monotonicity. At the junctions of three or more branches, the process advances with the minimum geodesic time while maintaining the coherence of the zero-level set. To constrain the computational domain, a narrow band of fixed thickness is defined:

[0133]

[0134] Where w is the narrowband thickness, a constant given by the voxel physical dimensions.

[0135] To suppress the numerical drift that accumulates over time after initialization, for f i (x) Periodically reinitialize according to the standard reset equation, preserving the properties of the distance function and not shifting the zero level set:

[0136]

[0137] φ(x,0)=f i (x)

[0138] Where φ is the temporary range field, τ is the pseudo-time, and sign(·) is the sign function. Reinitialization only occurs when... The process is performed internally, and f is replaced with φ after convergence. i .

[0139] Level set updates are performed within the narrowband, driving the evolution of the zero level set with normal velocity:

[0140]

[0141] Where t is the evolution time, δ τ (·) is the smooth Dirac function, v i (x) represents the normal velocity. The normal velocity consists of the data alignment term and the shape derivative of the geometric regularity, and can be superimposed with the bias voltage formed by the mass-conserved dual variables and the non-overlap penalty driven by the adjacency relationship. This invention only undertakes the general update form; the specific velocity composition has been given in adjacent chapters and will not be repeated. The time step satisfies the Courant condition to ensure stability, and the spatial derivative adopts the upwind scheme to avoid tangential numerical oscillations. After every few update steps, f is... i (x) Call a reinitialization once, then refresh synchronously. Geometric quantities such as normal and curvature.

[0142] To maintain compatibility among multiple cells, a consistent zero-level set signature is used at the common boundary voxels during initialization; if different units have different f... i (x) When a conflict occurs in the outer domain, the distance sign corresponding to the power graph affiliation is preferentially retained, and the conflict is automatically resolved by the energy term in the subsequent non-overlapping penalty. For extremely small isolated segments, the nearest neighboring unit is incorporated using a distance threshold strategy without changing the zero-level set, in order to reduce topological noise in subsequent evolution.

[0143] The benefits of this step include: first, the signed distance function strictly conforms to the anisotropic geometry of the spatial metric tensor field, and subsequent normal advance and curvature calculations are consistent with changes in imaging resolution; second, the fixed-thickness narrow band restricts the computation to the neighborhood, significantly reducing computational load and improving numerical stability; and third, periodic reinitialization preserves... Approaching 1 avoids distance degradation after long-range evolution, thus ensuring the reliability of the subsequent energy descent direction.

[0144] Example: On a set of samples containing attached structures, the initial boundary is first defined using a set of three-dimensional unit voxels. The anisotropic Eikonal equation is then solved using a spatial metric tensor field to obtain f. i (x). Set the narrowband thickness to a certain number of voxels, and construct... Within the narrow band, upwind schemes and Courant constraints are used for level set updates. The normal velocity is composed of a combination of the previously obtained data alignment term and geometric regularization. Reinitialization is performed every 10 to 20 time steps. During evolution, the zero level set bends along the high-cost direction of the microcarrier normal in the attached region without crossing the solid interface. In the deep region, a smooth transition band consistent with the resolution is obtained due to the anisotropy metric. If local isolated small fragments appear, they are incorporated into neighboring cells at the narrow band boundary using a distance threshold. After iterative convergence, f satisfies the distance function property. i (x) Stable zero-level sets and voxel-level masks provide a reliable geometric starting point for subsequent joint optimization and closed-loop verification.

[0145] Preferably, the data alignment term includes a boundary alignment term constructed based on the cell membrane probability map and intensity gradient direction when cell membrane channels are present, and a region consistency term constructed based on the intracellular intensity model and the extracellular intensity model when only cytoplasmic channels are present. The data alignment term drives the evolution of the signed distance function with normal velocity within the narrow band.

[0146] This invention establishes a differentiable registration relationship between observed intensity and geometric boundary within a narrow band of fixed thickness, using the normal velocity of a signed distance function as the carrier. Depending on the available channels, a boundary alignment term or a region consistency term is constructed; these can be enabled individually or in a weighted combination. Ultimately, the normal velocity drives the evolution of the signed distance function, outputting a boundary position consistent with the channel information.

[0147] When cell membrane channels exist, a boundary alignment term based on the cell membrane probability map and the intensity gradient direction is used. The cell membrane probability map is defined as P. mem (x), Where x is the spatial location, I mem (x) represents the cell membrane channel strength, g mem (x) is its spatial gradient. The unit normal of the signed distance function is: The normal velocity induced by the boundary alignment term is then written as:

[0148]

[0149] Where, α mem With β memLet be the weights, <·,·> be the inner product, and sign(·) be the sign function. This velocity term ensures that the zero-level set is attracted synchronously along the direction of membrane probability increase and the direction of intensity transition normal, avoiding tangential drift. To reduce the misleading effect of low-confidence regions, an uncertainty-based suppression weight is introduced:

[0150]

[0151] Where U(x) is the previously defined uncertainty field. The final boundary alignment velocity used for evolution is w. snr (x)v edge,i (x).

[0152] When only cytoplasmic channels exist, a region consistency term based on an internal-external statistical model is used. Let the cytoplasmic channel strength be denoted as ρ. cyt (x), the inner and outer log-likelihood costs are:

[0153]

[0154] Where, μ in With μ out σ is the internal and external mean. in With σ out The standard deviations are internal and external. The shape derivative of the regional consistency term gives the normal velocity at the zero level set. This means that when the cost on the outside of the boundary is greater than the cost on the inside, the boundary pushes outward; otherwise, it contracts inward. Statistical parameters are updated with a weighted average of the inside and outside voxels of the current geometry.

[0155]

[0156] Among them, Ω i For the cellular domain, The complement is then provided. To suppress the influence of low confidence, the above integral can be multiplied by w. snr (x) is used as the weight.

[0157] When cell membrane channels and cytoplasmic channels coexist, a weighted combination of observations is used to drive velocity:

[0158] v data,i (x)=λ edge w snr (x)v edge,i (x)+λ reg v reg,i (x)

[0159] Where, λ edge With λ reg Non-negative weights are assigned based on channel quality and batch consistency. Data alignment is calculated only within a fixed-thickness narrow band and is determined through... Participating in the normal progression of the zero-level set, where δτ (·) represents the smoothed Dirac function, with the same meaning as in the previous statement.

[0160] The following key points should be noted in the implementation of this invention. First, gradient calculation employs upwind difference and anisotropic filtering to ensure consistency with the spatial metric tensor field. Second, statistical parameter updates and velocity evaluation are performed alternately, with parameters re-estimated after several evolution steps each time to avoid jitter. Third, the velocity term is preferentially evaluated on the common boundary voxel set given by adjacency relationships, reducing computation in irrelevant regions. Fourth, regions with strong reflection or saturation are automatically weighted with uncertainty weights to prevent spurious boundary manipulation.

[0161] The effects of the above design are reflected in three aspects. First, the boundary alignment term ensures that the zero-level set is stably adsorbed at a position consistent with the probability and intensity transitions of the cell membrane, avoiding tangential drift and excessive jaggedness. Second, the region consistency term can still stably locate the cytoplasmic boundary through statistical differences between the inside and outside of the membrane even in the absence of membrane signals, maintaining a consistent shape at low signal-to-noise ratio levels. Third, uncertainty suppression naturally weakens both types of drivers in low-confidence regions, working in conjunction with the preceding spatial metric tensor field to avoid crossing microcarriers or traversing shadow regions.

[0162] Example: On a set of in vivo data containing membrane staining and cytoplasmic staining, the signed distance function and narrow band are first initialized using the previous method. The cell membrane probability map and intensity gradient are calculated to obtain the boundary alignment velocity; simultaneously, the mean and standard deviation of the inside and outside voxels of the current geometry are estimated to obtain the region consistency velocity. Based on the channel signal-to-noise ratio, a combination weight is set, and after superimposing uncertainty suppression, an observation-driven velocity is formed. Within the narrow band, the zero-level set is advanced at a Coulomb-stable time step, and the distance function is reinitialized and the statistical parameters are refreshed every 15 steps. After about several dozen rounds, the boundary stabilizes near the high cell membrane probability band and intensity gradient peak, the position with significant differences between the inside and outside of the main cytoplasmic region is consolidated, the adherent region does not cross the microcarrier, the high-depth layer maintains a smooth shape, and the output voxel-level mask is consistent with the subsequent closed-loop verification.

[0163] Preferably, the non-overlapping constraints include implementing a voxel-level smoothing penalty term at the common boundary of adjacent cells and performing an augmented Lagrangian update in combination with non-overlapping dual variables, and applying multi-body repulsion based on hypergraph relations at locations with tri- or multi-way contacts; the topological constraints include prohibiting the generation of holes and maintaining single connectivity during evolution.

[0164] This invention applies both non-overlapping and topological constraints within a narrow band, ensuring that multicellular boundaries remain independent, simply connected, and free of holes during evolution. Non-overlapping constraints, with the common boundary voxels as their domain, suppress geometric overlap through augmented Lagrangian forms; topological constraints, based on Eulerian characteristic and connectivity conservation, prevent the formation of holes and breaks. Both types of constraints only alter the normal velocity and duality, without changing the definitions of preceding metrics and data alignment terms.

[0165] Non-overlapping constraints are measured by overlap measure O ij =∫H(-f i (x))H(-f j (x))dx is the core, where x is the spatial location, f i (x) is the i-th signed distance function, and H(·) is the smoothed Heaviside function. Augmented Lagrange energy is written as follows:

[0166]

[0167] Where, μ ij For non-overlapping dual variables, γ aug This is the augmentation coefficient. For f i The shape derivative of (x) gives the non-overlapping normal velocity term:

[0168]

[0169] in, Let μ be the set of indices adjacent to the i-th cell. The dual variable μ is updated on a rolling basis according to the augmented Lagrange rule. ij ←μ ij +γ aug O ij This ensures that overlapping voxels are continuously penalized and gradually cleared during iteration.

[0170] At locations where ternary or multi-ary contacts exist, a hypergraph relation of many-body repulsion is introduced. The ternary contact measure is defined as:

[0171]

[0172] Multibody Repulsion Energy Writing in, Let η be the set of ternary groups extracted from adjacency relations. multi For weights. The corresponding additional term for normal velocity is:

[0173]

[0174] This causes the ternary junctions to contract into paired contacts, avoiding local "over-occupancy".

[0175] Topological constraints are based on the principles of simple connectivity and absence of holes. The "simple voxel" criterion of digital topology is used for candidate update points in each normal advance. Let the candidate position be x. * If f i (x *If changing from positive to negative or vice versa alters the number of connected components in the inner or outer domain, or changes the Euler characteristic χ, then the update is rejected; otherwise, the update is allowed. In implementation, a disjoint-set data structure with inner and outer labels is maintained within a narrowband, and "simple voxels" are determined in constant time using a lookup table for a local 3×3×3 neighborhood. When consecutively rejected topological events occur, the weights of curvature regularization and non-overlap penalties are temporarily increased to alleviate geometric conflicts and guide the boundary away from high-risk configurations. Synthetic normal velocity writing:

[0176]

[0177] Among them, v data,i (x) represents the data alignment speed, α dist With α curv κ represents the weight. i (x) is the mean curvature, v mass,i (x) represents the mass bias velocity. The evolution equation remains:

[0178]

[0179] Where, δ τ (·) represents the smoothed Dirac function. Non-overlapping and many-void repulsion are evaluated only on the set of adjacent common boundary voxels, avoiding globally independent computations; the augmentation coefficient γ... aug A gradual ascent strategy is adopted to achieve smooth convergence.

[0180] This design yields three direct benefits. First, the non-overlapping constraint transforms geometric conflicts into differentiable velocity corrections and dual updates, rapidly eliminating overlapping voxels without disrupting aligned boundaries. Second, multi-voxel repulsion suppresses intersections of ternary or higher units under complex contact topologies, reducing the frequency of subsequent splitting and merging. Third, topological constraints ensure that each cell domain remains simply connected and free of holes, making statistics and subsequent counterfactual tests stable and reliable.

[0181] Example: On samples containing high-density adherent regions, a signed distance function is first initialized using a power-graph 3D unit, and data-driven evolution is performed within a narrow band. After detecting several triangular contacts, a hypergraph ternary clique set is constructed, and multi-body repulsion is enabled; simultaneously, overlap measures are evaluated on the common boundary voxel set, and dual variables are updated in an augmented Lagrangian form. After dozens of iterations, the overlapping voxel count drops to zero, triangular contacts shrink into paired contacts, and any candidate topological changes are intercepted by the "simple voxel" criterion. The final output cell domain remains simply connected, pore-free, and pairwise non-overlapping, with consistent boundary and channel information and spatial metric tensor fields, allowing direct entry into the imaging consistency closed-loop and statistical stages.

[0182] Based on the synthetic density generated by the cell domain, synthetic observations are obtained through imaging operators and residuals are calculated. Based on the residuals, splitting, merging, or deformation is triggered and the power graph weights and spatial metric tensor field are updated. The process is iterated until convergence, and the final voxel-level mask is output.

[0183] Preferably, the synthesis density generated based on the cell domain includes constructing a thin-shell profile for the cell membrane and an internal profile for the cytoplasm on the normal coordinates of the cell domain. The residual between the synthetic observation and the original volume data is mapped to the correction amount of the boundary normal velocity and the correction amount of the power graph weight through the accompanying imaging operator, and splitting, merging or deformation is performed accordingly until convergence.

[0184] This invention constructs a closed loop for consistency between geometry and imaging physics: A synthetic density is generated based on the obtained cellular domain; synthetic observations are obtained through imaging operators and compared with volume data to form residuals; the residuals are fed back as boundary normal velocity corrections and power graph weight corrections using the adjoint function of the imaging operators; simultaneously, splitting, merging, or deformation are triggered according to the spatiotemporal structure of the residuals until convergence, outputting the final voxel-level mask. This closed loop maintains the same coordinates and metrics as the preceding continuous fluorescence density field, spatial metric tensor field, three-dimensional cells, and signed distance function, avoiding cross-domain inconsistencies.

[0185] To generate a synthesis density consistent with the channel within the cellular domain, normal coordinates of the zero-level set are used. Normal coordinates are defined as a signed distance function, and are written as follows for the cell membrane profile and cytoplasm profile:

[0186]

[0187] Where x is the spatial position, f i (x) is the signed distance function of the i-th cell, k shell For thin-shell kernel functions, σ i k is the thickness parameter of the thin shell. in For the internal profile kernel function, τ i For the internal attenuation scale, a i With c i Here is the amplitude parameter. The total synthesis density is...

[0188] Let the imaging operator and its companion be denoted as Where y is volume data, y syn For composite observations, r is the residual and g is the associated propagation field. Forward imaging operator, As the adjoint operator, θ is the set of physical parameters. The propagation field is transformed into a normal velocity correction at the zero level set:

[0189]

[0190] fi (x)=0

[0191] Where, γ v This is the weight. This term, together with the data alignment term, curvature regularization, mass bias, and non-overlap penalty, constitutes the synthetic normal velocity, achieving physical consistency correction of the geometric boundary.

[0192] To feed back the bias in voxel assignment from the observation residuals to the mass conservation layer, a weight correction is defined:

[0193]

[0194] Where, Δw i γ is the power graph weight adjustment amount. w For weights, L i (w) represents the i-th 3D unit of the current power graph. The triggering of this correction, splitting, merging, and deformation is based on the spatial structure and adjacency of the residuals. For each cell domain, if a stable bimodal residual exists within the domain and is separated by a saddle-shaped valley, and the mass conservation dual variable is consistently positive while the residual outside the boundary is negative, then a candidate splitting surface is proposed at the saddle-shaped valley, and a new zero-level set is initialized with a small perturbation; if two adjacent domains have complementary residual signs on both sides of the common boundary, and both the boundary curvature and mass bias point towards merging, then the common boundary is revoked and merging is performed; when only a unilateral residual or slight misalignment occurs, deformation takes priority, through Δv i The local advancement of (x) completes the correction. All the above topology operations comply with the aforementioned topology constraints and non-overlapping constraints to avoid introducing holes or multiple overlaps.

[0195] The spatial metric tensor field is updated using residual-driven uncertainty modulation. Let the uncertainty field be U(x), then it is updated according to a monotonically increasing function ψ:

[0196] U(x)←ψ(|r(x)|)

[0197]

[0198] Where G(x) is the spatial metric tensor field, G psf (x) represents the anisotropic component obtained from the point spread function, τ is the uncertainty weight, d(x) and n(x) are the distance to the microcarrier surface and the normal vector, respectively, κ is the normal cost gain, and I is the unit tensor. The metric is updated only on the uncertainty component and does not change the physical components of the point spread function and depth decay.

[0199] Convergence is determined using multiple criteria: the residual norm decreases to a threshold and its change is less than a set proportion over several consecutive rounds; the quality error of all cells meets the tolerance and the quality-conserving dual variable no longer drifts significantly; the zero-level set has a maximum normal displacement less than a threshold and no new division or merging events. Once these conditions are met, the cell domain is solidified and a voxel-level mask is output, along with statistical indicators such as cell number, volume, and adherent contact area.

[0200] Example: On samples containing both adherent and suspended cells, a thin shell and internal profile are first constructed for each cell according to normal coordinates to form a synthetic density, which is then used by an imaging operator to generate synthetic observations. Residuals are calculated, and the return field is obtained from the accompanying return. A normal velocity correction is formed at the zero level set, which, along with data alignment terms, curvature regularization, mass bias, non-overlap, and multi-body repulsion, advances the boundary. Simultaneously, the power graph weights are updated using the voxel integral of the return field, and the uncertainty components are updated using the residual amplitude. After several dozen rounds, the boundary of the adherent region bends along the high-cost direction of the microcarrier normal but does not cross it; the deep weak signal region obtains a smooth boundary through mass and accompanying corrections; some excessively large cells are split under the cues of residual bimodalities and dual variables, while adjacent excessively fine cells merge under the cues of complementary residuals. Finally, the residuals stabilize, mass conservation is satisfied, and the boundary displacement tends to zero. The output voxel-level mask is consistent with the channel observations, providing repeatable quantitative results for batch comparison and quality release.

[0201] The above are merely embodiments of this application and are not intended to limit the scope of this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of the claims of this application.

Claims

1. A method for three-dimensional segmentation of umbilical cord stem cells, characterized in that, Includes the following steps: A continuous fluorescence density field is generated by inverting the confocal imaging or light sheet imaging volume data according to the imaging operator. The continuous fluorescence density field is then used to generate the initial synthetic observation in the forward direction through the imaging operator. The initial residual is obtained based on the difference between the initial synthetic observation and the confocal imaging or light sheet imaging volume data. A spatial metric tensor field is constructed based on the point spread function, depth attenuation, and uncertainty terms generated from the initial residual and signal-to-noise ratio. The spatial metric tensor field is a spatial position-related tensor field used to characterize the cost of motion in different directions in three-dimensional space. Under the geodesic distance induced by the spatial metric tensor field, the cell center and target mass are determined based on the continuous fluorescence density field. A power graph is constructed by dividing the data using a power distance function. The power distance function is obtained by subtracting the power graph weight corresponding to each cell center from the square of the geodesic distance. The power graph weight is used to adjust the belonging boundary of the corresponding three-dimensional unit. The power graph weight is optimized by semi-discrete optimal transport to obtain the three-dimensional unit and adjacency relationship that satisfy the mass conservation. The signed distance function is initialized with three-dimensional units under a spatial metric tensor field. Data alignment terms are constructed based on the continuous fluorescence density field and the signed distance function is jointly optimized. Non-overlapping constraints and topological constraints are applied to obtain cell domains and voxel-level masks. Based on the synthetic density generated by the cell domain, the segmented synthetic observation is obtained through the imaging operator and the segmentation residual is calculated. Based on the segmentation residual, splitting, merging or deformation is triggered, and the power graph weights and the uncertainty term in the spatial metric tensor field are updated according to the segmentation residual. The process is iterated until convergence, and the final voxel-level mask is output.

2. The method according to claim 1, characterized in that, Constructing the spatial metric tensor field involves forming an anisotropy metric based on the principal direction and principal axis length of the point spread function, determining the resolution variation at different depths based on depth attenuation, generating an uncertainty term based on the initial residual and the signal-to-noise ratio estimated by channel, and combining the anisotropy metric, uncertainty term, and unit tensor to obtain the spatial metric tensor field. In the presence of microcarriers, the motion cost in the normal direction is increased at positions close to the microcarrier surface based on the microcarrier surface distance field and microcarrier normal to form an impassable constraint. Inverting and generating the continuous fluorescence density field involves performing flat-field correction, fringe suppression, and depth attenuation normalization on the volume data, and solving the non-negative regularized inverse problem based on the imaging operator.

3. The method according to claim 1, characterized in that, Determining the cell center involves constructing a superlevel set sequence on the nuclear channels of a continuous fluorescence density field, screening for local maxima as candidate points based on topological persistent homology, and performing nonmaximum suppression and nearest neighbor merging under geodesic distance induced by a spatial metric tensor field to obtain the cell center.

4. The method according to claim 1, characterized in that, Determining the target quality involves volume integration of cytoplasmic channels in a continuous fluorescence density field within a geodesic neighborhood induced by a spatial metric tensor field with the cell center as the nucleus, and amplitude limiting based on the value range determined by historical batch statistics. In the absence of cytoplasmic channels, the target quality is determined by the intensity integration of available channels based on the calibrated channel mapping.

5. The method according to claim 1, characterized in that, Constructing a power graph involves: solving the geodesic distance field from each cell center to the voxel position under the spatial metric tensor field; obtaining the power distance function by subtracting the power graph weight of the corresponding cell center from the square of the geodesic distance; determining the voxel affiliation according to the cell center corresponding to the minimum power distance; generating the power graph and obtaining the adjacency relationship.

6. The method according to claim 1, characterized in that, The optimization of semi-discrete optimal transport involves convex optimization iteration of the power graph weights, making the volume fraction of the continuous fluorescence density field in the cytoplasmic channel equal to the corresponding target mass, and generating a mass-conserving dual variable for the mass bias of the signed distance function.

7. The method according to claim 1, characterized in that, Initializing the signed distance function involves solving the Eikonal equation for distance calculation at the boundary of the 3D element to obtain the initial signed distance function, and updating the level set within a narrow band of fixed thickness.

8. The method according to claim 1, characterized in that, The data alignment term includes: when cell membrane channels exist, a cell membrane probability map is generated based on the cell membrane channels. The cell membrane probability map is used to represent the confidence distribution of each voxel belonging to the cell membrane boundary, and a boundary alignment term is constructed in combination with the intensity gradient direction of the cell membrane channels; when only cytoplasmic channels exist, a region consistency term is constructed based on the intracellular intensity model and the extracellular intensity model; the data alignment term drives the evolution of the signed distance function with normal velocity within the narrow band.

9. The method according to claim 1, characterized in that, Non-overlapping constraints include: applying voxel-level smoothing penalties to the common boundaries of adjacent cells and performing augmented Lagrangian updates in conjunction with non-overlapping dual variables; constructing hypergraph relations and applying multi-body repulsion based on hypergraph relations when there are ternary or multi-ary contacts, wherein the hypergraph relations are formed by cell domains as nodes and by three or more cell domains involved in the same ternary or multi-ary contact location; topological constraints include prohibiting the generation of holes and maintaining single connectivity during evolution.

10. The method according to claim 1, characterized in that, Based on the cell domain, the generated synthetic density includes constructing a thin-shell profile of the cell membrane and an internal profile of the cytoplasm on the normal coordinates of the cell domain. The segmentation residual between the segmented synthetic observation and the confocal imaging or light sheet imaging volume data is mapped to the correction amount of the boundary normal velocity and the correction amount of the power graph weight through the accompanying imaging operator. Based on this, splitting, merging or deformation is performed until convergence.

Images