Organoid automatic labeling method based on machine vision and deep learning

Abstract

The application relates to the technical field of image recognition processing, in particular to an organoid automatic labeling method based on machine vision and deep learning, which comprises the following steps: acquiring an image to be processed under a high-throughput culture environment and performing image preprocessing; inputting the preprocessed image to be processed into a preset digital image processing framework; the digital image processing framework comprises an image filtering unit, a three-dimensional boundary reconstruction unit and a texture gray scale analysis unit; an independent target region mask is generated through the three-dimensional boundary reconstruction unit; a high-density core mask is generated through the texture gray scale analysis unit; a three-dimensional geometric topology graph of the organoid is constructed based on the target region mask and the high-density core mask; and a three-dimensional space form feature of a neighborhood node is obtained through processing by using a numerical iterative aggregation algorithm; and the three-dimensional space form feature is analyzed to output an image feature parameter.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of image recognition processing, in particular to an organoid automatic labeling method based on machine vision and deep learning. BACKGROUND

[0002] In the field of biomedical engineering and drug screening, organoids, as a kind of three-dimensional cell culture model constructed in vitro, can simulate the structure and function of real organs. In order to evaluate the drug response or growth state, it is usually necessary to use high-throughput microscopic imaging system to collect a large amount of organoid bright field images. However, such images have very high data complexity, which brings significant challenges to traditional digital image processing technology.

[0003] Firstly, in terms of image segmentation and boundary extraction, the high-throughput culture environment leads to a very high distribution density of target objects in the image. In two-dimensional projection images, multiple three-dimensional targets often have physical contact, overlap or occlusion, forming complex connected domains. Traditional methods based on gray threshold, edge detection or conventional watershed algorithm are difficult to accurately separate the dense adhesion area from the pixel level according to the weak gradient change, often leading to under-segmentation or contour distortion of the target in the image processing result.

[0004] Secondly, in terms of texture feature analysis, the target object often contains complex heterogeneous structure, such as high-density necrotic core. The gray difference between these areas and the surrounding healthy tissue fluctuates greatly in different imaging batches, and is easily disturbed by imaging noise. Existing pixel statistical methods lack spatial topological constraints, and are difficult to accurately extract specific texture gray areas inside the target while suppressing background noise.

[0005] In addition, existing image analysis schemes are mostly limited to morphological measurement of discrete two-dimensional targets, such as calculating area and perimeter, ignoring the three-dimensional geometric topological relationship of target objects in the real physical space. Two-dimensional image data loses the spatial adjacency and volume compression information between targets, making it difficult to restore the group distribution state and microenvironment interaction characteristics of target objects from image data.

[0006] Therefore, there is an urgent need for an image processing method that can integrate spatial boundary reconstruction, texture statistical analysis and geometric topological calculation to solve the problems of adhesion segmentation difficulty, internal feature extraction inaccuracy and spatial information loss in high-density microscopic images.

[0007] Therefore, an organoid automatic labeling method based on machine vision and deep learning is proposed. SUMMARY

[0008] The application aims to provide an organoid automatic labeling method based on machine vision and deep learning, which generates independent target region masks and high-density core masks through the digital image processing architecture, constructs a three-dimensional geometric topology graph of the organoid, analyzes the three-dimensional spatial morphological features, outputs image feature parameters, and realizes accurate image labeling.

[0009] To achieve the above-mentioned purpose, the application provides the following technical scheme.

[0010] The organoid automatic labeling method based on machine vision and deep learning comprises the following steps:

[0011] An image to be processed in a high-throughput culture environment is obtained, and image preprocessing is performed.

[0012] The preprocessed image to be processed is input into a preset digital image processing architecture; the digital image processing architecture comprises an image filtering unit, a three-dimensional boundary reconstruction unit, and a texture gray scale analysis unit.

[0013] Through the three-dimensional boundary reconstruction unit, pixel clustering based on a gradient vector field is performed, a centripetal unit vector field of each pixel point in the image pointing to the local geometric center to which the pixel point belongs is calculated, and dense adhesion between connected domains is numerically segmented based on the divergence characteristics of the vector field, thereby generating independent target region masks.

[0014] Through the texture gray scale analysis unit, a pixel gray scale statistical analysis task of the image is performed; under the constraint of topology inclusion, the gray scale abnormal region of the target region is extracted, and a high-density core mask is generated.

[0015] Based on the target region mask and the high-density core mask, a three-dimensional geometric topology graph of the organoid is constructed; the three-dimensional spatial morphological features of the neighborhood nodes are obtained by using a numerical iterative aggregation algorithm, the three-dimensional spatial morphological features are analyzed, and image feature parameters are output.

[0016] The design of the digital image processing architecture comprises:

[0017] The image filtering unit fuses low-frequency global contour features and high-frequency local detail features, adapts to the scale change of the target region at different stages, and outputs a multi-scale filtering response map.

[0018] The three-dimensional boundary reconstruction unit processes gradient edges and low-frequency contour information in the image, and constructs geometric shapes.

[0019] The texture gray scale analysis unit processes high-frequency texture details and gray scale statistical features in the image, calculates and separates the gray scale differences of different pixel density regions.

[0020] The target region mask acquisition process comprises:

[0021] A normalized direction vector field is constructed by calculating the normalized direction vector of each pixel in the image pointing to the local geometric center to which the pixel belongs; meanwhile, a membership value of each pixel belonging to the target foreground region is calculated to construct a foreground probability distribution map;

[0022] A vector field divergence map is generated by performing divergence operation on the normalized direction vector field; the divergence operator is used to detect the mutation area of the vector direction field, and the positive peak area of the divergence is defined as the adhesion boundary barrier potential between target regions;

[0023] Hough voting is performed on the normalized direction vector field, the convergence center of the vector flow is counted, and the area with convergence intensity exceeding a preset threshold is determined as a geometric center marker, which is used as a water injection seed for the watershed algorithm;

[0024] The terrain potential surface is constructed by superimposing the foreground probability distribution map and the vector field divergence map with the geometric center marker as the starting point, and a water injection flooding operation is performed to generate a mask of independent target regions.

[0025] The terrain potential surface acquisition process includes:

[0026] The foreground probability distribution map is subjected to a numerical inversion operation to construct a basic terrain map, so that the internal area of the target region is mapped as a low potential basin, and the background region is mapped as a high potential highland;

[0027] The vector field divergence map is subjected to linear rectification processing to filter out negative values and retain positive areas, extract adhesion boundary features representing the vector departure trend, and construct a boundary barrier potential map;

[0028] The basic terrain map and the boundary barrier potential map are added at the pixel level to generate the final terrain potential surface.

[0029] The process of identifying the gray abnormal area by the texture gray analysis unit includes:

[0030] The multi-scale filtering response map output by the image filtering unit is received by the texture gray analysis unit; a multi-scale feature fusion operation is performed to splice the features containing high-resolution details in the shallow layer and the features containing semantic information in the deep layer, and a feature response sensitive to local optical density changes of the image is extracted through a spatial filtering matrix to obtain a texture difference feature representing the difference between the gray aggregation area inside the target region and the surrounding gray smooth area;

[0031] Based on the extracted texture difference feature, a pixel-level numerical mapping is performed to output a high-density area distribution map, wherein the pixel value represents the confidence degree of belonging to the gray abnormal area;

[0032] The target region mask is used as a spatial constraint condition to perform a spatial intersection operation with the high-density region distribution map to filter out noise pixels located outside the individual contour of the target region, and to ensure that the high-density core is contained within the target region in terms of spatial topology.

[0033] The high-density region distribution map after topological filtering is subjected to Gaussian smoothing processing to fit the fuzzy boundary characteristics of the abnormal gray region, and an adaptive threshold is set for binary processing to generate a high-density core mask.

[0034] The process of constructing a three-dimensional geometric topology map includes:

[0035] The generated target region mask is subjected to connected domain labeling, and each independent connected domain is instantiated as a grid node to construct a set of grid nodes.

[0036] Based on the target region mask, the geometric morphological characteristics of each grid node are calculated, including at least the equivalent diameter and the circularity; and based on the high-density core mask, the pixel statistical characteristics of each node are calculated, including at least the area ratio of the abnormal gray region and the average gray level of the abnormal gray region; the geometric morphological characteristics and the pixel statistical characteristics are combined to generate an initial feature vector of each node.

[0037] The Euclidean distance between the centroids of the target regions corresponding to any two grid nodes is calculated; if the Euclidean distance is less than a predetermined interaction threshold, it is determined that there is spatial interaction between the two, and a connection edge is established between the corresponding two nodes, thereby generating a three-dimensional geometric topology map including a node set, an edge set, and a node feature vector.

[0038] The process of processing the three-dimensional geometric topology map using a numerical iterative aggregation algorithm includes:

[0039] The three-dimensional geometric topology map is input into a three-dimensional topology calculation module, and for each target node, the data similarity between the target node and all connected neighbor nodes is calculated to generate a normalized correlation coefficient representing the influence weight of the neighbor nodes on the three-dimensional spatial state of the target node.

[0040] Based on the correlation coefficient, the initial feature vector of the neighbor node is subjected to a weighted summation operation to transfer the morphological crowding information and the gray distribution information of the neighbor to the target node to generate an aggregated feature that integrates environmental context information.

[0041] The aggregated feature is input into a spatial mapping matrix for feature transformation to update the feature vector of the target node and obtain a three-dimensional spatial morphological feature.

[0042] The process of analyzing the three-dimensional spatial morphological feature includes:

[0043] All nodes in the three-dimensional geometric topology graph are traversed, and the local physical crowding degree of each target region is calculated based on the number of connection edges of each node in combination with the distance weight of the neighborhood node;

[0044] The updated node feature vector, local physical crowding degree and area proportion of the gray abnormal region are taken as joint input features; through feature space clustering analysis, the distribution difference of different gray modes in the feature space is mined to judge the state type of the target region; and the image feature parameters are obtained through the geometric morphological features of the grid node and the state type.

[0045] Compared with the prior art, the beneficial effects of the present application are:

[0046] 1. The decoupling design of "shared backbone and independent branch" is adopted in the present application, and the multi-scale and multi-task conflict problem in organoid image analysis is effectively solved. Through the feature pyramid structure, the edges of small cell groups can be accurately outlined by using shallow features, and the overall morphology of large organoids can be accurately recognized by using deep features, so that missing detection is avoided. Boundary reconstruction focuses on the low-frequency contour and gradient change of the image, and texture analysis focuses on the high-frequency details and gray statistics of the image; the two are separated for learning in different branches, avoiding feature interference when a single network processes a target with clear contour but complex internal structure, and significantly improving the segmentation accuracy in a complex background.

[0047] 2. The field theory idea in physics is creatively introduced to solve the image segmentation problem in the present application; the traditional segmentation method based on binary mask can only be hard cut by morphological erosion when processing adherent targets, which can easily lead to target shape distortion or segmentation failure. The topographic potential surface constructed by the present scheme skillfully utilizes the sudden change characteristics of the vector direction at the adherent position; even if two organoids have almost no difference in gray scale, as long as their geometric centers are different, the vector field will inevitably deviate at the contact surface. This segmentation method based on mathematical potential surface has high robustness for high-density and overlapping organoids, and the segmentation accuracy is significantly better than that of the traditional method.

[0048] 3. The present application breaks through the limitation of two-dimensional image analysis, and restores the physical occupation relationship of organoids in the real culture environment through the "sphere assumption". The traditional two-dimensional distance calculation often ignores the influence of object size, and large objects are more likely to contact than small objects. The constructed topology graph truly reflects the occupation relationship, providing a key structured data basis for subsequent differentiation of necrosis caused by physical crowding. BRIEF DESCRIPTION OF DRAWINGS

[0049] Figure 1 The flowchart of the organoid automatic labeling method based on machine vision and deep learning of the present application;

[0050] Figure 2A data logic diagram of the organoid automatic labeling method based on machine vision and deep learning of the present application is shown in the figure;

[0051] Figure 3 A structural diagram of the digital image processing architecture of the present application is shown in the figure. DETAILED DESCRIPTION

[0052] The technical solutions in the embodiments of the present application will be clearly and completely described below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative work fall within the scope of protection of the present application.

[0053] Embodiment one:

[0054] The present application proposes an organoid automatic labeling method based on machine vision and deep learning. The flow of the method is shown in the figure Figure 1 The data logic is shown in the figure Figure 2 , which includes:

[0055] The image to be processed under the high-throughput culture environment is acquired, and image preprocessing is performed.

[0056] The preprocessed image to be processed is input into a preset digital image processing architecture. The digital image processing architecture includes an image filtering unit, a three-dimensional boundary reconstruction unit, and a texture gray scale analysis unit.

[0057] Through the three-dimensional boundary reconstruction unit, pixel clustering based on the gradient vector field is performed, the centripetal unit vector field of each pixel point in the image pointing to the local geometric center is calculated, and the dense adhesion between connected domains is numerically segmented based on the divergence characteristics of the vector field, to generate an independent target region mask.

[0058] Through the texture gray scale analysis unit, the pixel gray scale statistical analysis task of the image is performed. The gray scale abnormal area of the target region is extracted under the topological inclusion constraint, and a high-density core mask is generated.

[0059] Based on the target region mask and the high-density core mask, a three-dimensional geometric topology graph of the organoid is constructed. The three-dimensional spatial morphological characteristics of the neighborhood nodes are obtained by using a numerical iterative aggregation algorithm, and the three-dimensional spatial morphological characteristics are analyzed to output the image feature parameters.

[0060] Preferably, the image acquisition and preprocessing process includes:

[0061] An automatic high-throughput scanning microscope is used to acquire bright field images of organoids in a 96-well plate or a 384-well plate. Since there is often uneven illumination and background noise in microscopic imaging, the image is first preprocessed.

[0062] Specifically, the image is converted to a lightness-chrominance color space (such as Lab space), and the lightness channel is extracted; a limited contrast adaptive histogram equalization algorithm is applied to divide the image into multiple small blocks, and histogram equalization is performed on each block, and then bilinear interpolation is used to eliminate block effects; this step can significantly enhance the contrast of target edges and background, while suppressing high light noise.

[0063] Preferably, the structure of the digital image processing architecture is as shown in Figure 3 , comprising:

[0064] The image filtering unit fuses low-frequency global contour features and high-frequency local detail features, adapts to the scale changes of the target region at different stages, and outputs a multi-scale filtering response map.

[0065] The three-dimensional boundary reconstruction unit processes gradient edges and low-frequency contour information in the image to construct geometric shapes.

[0066] The texture gray scale analysis unit processes high-frequency texture details and gray scale statistical features in the image, and calculates and separates the gray scale differences of different pixel density regions.

[0067] It should be noted that the three-dimensional boundary reconstruction in the present application does not refer to the measurement of the physical thickness of the target object, but refers to the mathematical three-dimensional manifold structure constructed based on two-dimensional pixel coordinates and calculated potential eigenvalues in the digital image processing process. By performing watershed operation in the virtual three-dimensional space, numerical separation of the two-dimensional overlapping boundaries is achieved.

[0068] Specifically, constructing geometric shapes refers to constructing a mathematical three-dimensional terrain surface based on gradient potential, and segmenting geometric shapes based on the terrain surface.

[0069] The digital image processing architecture in the present embodiment adopts a backbone double-branch design:

[0070] Image filtering unit: a deep residual network (such as a 50-layer deep residual network) is used as the backbone.

[0071] In order to solve the scale difference of organoids from small cell clusters (several tens of microns) to mature bodies (several hundred microns), a feature pyramid structure is introduced. Specifically, the low-resolution, strong semantic feature map output by the deep layer of the backbone network is upsampled, and the high-resolution, strong detail feature map output by the shallow layer is connected horizontally and added element by element. Finally, four feature maps of different scales are output, corresponding to one-fourth, one-eighth, one-sixteenth, and one-thirty-second of the original image size, respectively.

[0072] Channel attention mechanism: before the feature map enters the subsequent processing unit, a channel attention module is embedded.

[0073] This module compresses spatial information through global average pooling, then learns the weight coefficients of each feature channel through two fully connected layers, and finally multiplies the original feature map back with the weight; so that the model can automatically pay attention to the feature channels more sensitive to boundaries or textures.

[0074] Three-dimensional boundary reconstruction unit: focus on geometric contour recovery; this unit contains convolutional layers and up-sampling layers, outputting a centripetal unit vector field and a foreground probability distribution map. This unit contains a series of convolutional layers and up-sampling layers, outputting two tensors:

[0075] Vector field prediction tensor: contains two channels, representing the horizontal and vertical displacement of the pixel pointing to the target center; foreground probability tensor: contains one channel, representing the confidence of the pixel belonging to the target foreground.

[0076] Texture gray scale analysis unit: focus on internal texture analysis; use dilated convolution to expand the receptive field, specifically respond to high-density texture areas. This unit uses dilated convolution technology, also known as dilated convolution; by inserting holes in the convolution kernel, the receptive field is expanded without increasing the number of parameters; this is to ensure that when detecting large areas of necrosis inside the organoid, the network can perceive a large enough context, rather than just focusing on local noise.

[0077] The training process of the digital image processing architecture includes:

[0078] Data preparation: collect micrograph images of organoids, and perform pixel-level labeling, including individual contours and internal necrotic areas.

[0079] Loss function: multi-task joint optimization is used for training.

[0080] For the vector field, use the smooth L1 loss function to calculate the distance error between the predicted vector and the true vector.

[0081] For the foreground probability and necrotic area, use the binary cross-entropy loss function combined with the Dice coefficient loss function to deal with the sample imbalance problem between foreground and background.

[0082] Training strategy: use the stochastic gradient descent optimizer to update parameters. Random rotation, mirror flip, and simulated out-of-focus Gaussian blur are used as data augmentation methods during training to improve the generalization ability of the model.

[0083] The application adopts a decoupling design of "shared backbone and independent branch", effectively solving the multi-scale and multi-task conflict problem in organoid image analysis. Through the feature pyramid structure, the shallow layer features can be used to accurately outline the edges of small cell groups, and the deep layer features can be used to accurately identify the overall morphology of large organoids, avoiding missing detection. Boundary reconstruction focuses on the low-frequency profile and gradient change of the image, while texture analysis focuses on the high-frequency details and gray scale statistics of the image; separating the two for learning in different branches avoids feature interference when a single network processes targets with clear outlines but complex interiors, significantly improving the segmentation accuracy in complex backgrounds.

[0084] Preferably, the target region mask acquisition process comprises:

[0085] A normalized direction vector of each pixel in the image pointing to the local geometric center to which it belongs is calculated to form a centripetal unit vector field. At the same time, the membership value of each pixel belonging to the target foreground region is calculated to form a foreground probability distribution map.

[0086] The centripetal unit vector field is subjected to divergence operation to generate a vector field divergence map. The divergence operator is used to detect the mutation area of the vector direction field, and the positive peak area of the divergence is defined as the adhesion boundary barrier potential between target regions.

[0087] The centripetal unit vector field is subjected to Hough voting, and the convergence center of the vector flow is counted. The region with a convergence intensity exceeding a preset threshold is determined as a geometric center marker, which is used as a water injection seed for the watershed algorithm.

[0088] The geometric center marker is used as a starting point, and the foreground probability distribution map is superimposed on the vector field divergence map to construct a terrain potential surface. A water injection flooding operation is performed to generate mutually independent target region masks.

[0089] Preferably, the terrain potential surface acquisition process comprises:

[0090] The foreground probability distribution map is subjected to a numerical inversion operation to construct a basic terrain map, so that the internal region of the target region is mapped to a low potential basin, and the background region is mapped to a high potential highland.

[0091] The vector field divergence map is subjected to linear rectification processing to filter out negative values and retain positive regions, extract adhesion boundary features representing the trend of vector departure, and construct a boundary barrier potential map.

[0092] The basic terrain map and the boundary barrier potential map are added at the pixel level to generate the final terrain potential surface.

[0093] Generation of centripetal unit vector field: for each pixel in the image, a two-dimensional vector pointing to its local geometric center is predicted by the regression prediction module. To unify the scale, the vector is normalized to have a length of 1, retaining only the directional information, thus forming a centripetal unit vector field. Inside the target, the vector flow lines converge; in the background area, the vector is zero or random.

[0094] Vector field divergence operation: calculate the divergence of the vector field, which is a scalar describing the flux of the vector field.

[0095] Calculation method: take the horizontal partial derivative of the horizontal component of the vector field, and take the vertical partial derivative of the vertical component. The sum of the two is the divergence. Inside the organoid, the vector points to the center, and the divergence is negative (convergent); on the boundary line between the two organoids, the left vector points to the left, and the right vector points to the right, the direction is away, and the divergence presents a large positive value (divergent).

[0096] Construction of mathematical three-dimensional terrain potential surface: in order to use the watershed algorithm, a virtual three-dimensional terrain needs to be constructed;

[0097] Construction of basin (basic terrain): perform a numerical negation operation on the foreground probability distribution map, for example, subtract the probability value by 1; the probability inside the target is close to 1, and after negation, it is close to zero, forming a low potential basin; the probability of the background is close to zero, and after negation, it is close to 1, forming a high potential highland;

[0098] Construction of ridge (barrier potential): linear rectification is performed on the divergence map to filter out negative values and only keep positive values, and a preset weight coefficient is multiplied to convert the divergent feature at the adhesion site into a towering "ridge";

[0099] Superposition: add the "basin" and "ridge" at the pixel level to generate the final terrain potential surface.

[0100] Seed point extraction: perform Hough voting on the vector field, and each pixel casts a vote along its vector direction. The peak value formed by the cumulative votes is the geometric center marker of the target.

[0101] Water filling: take the geometric center marker as the starting water injection point and perform water filling operation on the constructed terrain potential surface. The water flow will fill its own basin and stop when it encounters the "ridge" composed of divergence.

[0102] The specific Hough voting process is as follows: initialize an accumulator matrix consistent with the size of the original image. For each pixel in the predicted vector field, project within a preset radius, such as 5-50 pixels, along its vector direction; add 1 to the corresponding position of the accumulator. Finally, perform non-maximum suppression on the accumulator to extract the local peak points as the predicted geometric center coordinates.

[0103] The application creatively introduces the field theory idea in physics to solve the image segmentation problem; the traditional segmentation method based on binary mask can only be hard cut by morphological erosion when dealing with adherent targets, which is easy to cause target shape distortion or segmentation failure. The topographic potential surface constructed by the scheme skillfully uses the sudden change characteristics of the vector direction at the adherent position; even if two organoids have almost no difference in gray scale, as long as their geometric centers are different, the vector field will inevitably deviate at the contact surface. This segmentation method based on mathematical potential surface has high robustness for high-density and overlapping occluded organoids, and the segmentation accuracy is significantly better than that of traditional methods.

[0104] Further, the process of predicting by regression the normalized direction vector of each pixel in the image pointing to the local geometric center to which it belongs, specifically includes:

[0105] Constructing a regression prediction module: a pixel-level regression prediction module is further configured at the end of the three-dimensional boundary reconstruction unit, which outputs a prediction map containing two feature channels, wherein the first channel value represents the horizontal displacement component of the current pixel pointing to the geometric center of the target region, and the second channel value represents the vertical displacement component;

[0106] Performing numerical normalization processing: performing root mean square operation on the horizontal displacement component and the vertical displacement component to obtain the displacement module length, and dividing the displacement component by the displacement module length to generate a centripetal unit vector field with a module length of 1;

[0107] The local geometric center is determined in the model training stage by the following method: calculating the image moment of the manually labeled target region connected domain, and determining the geometric center coordinates of the connected domain by using the ratio of zeroth moment to first moment.

[0108] In order to accurately predict the vector field, first, a standard answer, i.e. the true value label, needs to be generated based on the manually labeled data:

[0109] Geometric center calculation: for each independently connected target region mask in the training set, the geometric center is calculated using the image moment algorithm. Specifically, the first moment of the target region in the horizontal and vertical directions is calculated, and then divided by the zeroth moment of the target region, i.e. the area; thus obtaining the unique geometric center coordinates of the target region, i.e. the geometric centroid.

[0110] Displacement vector generation: traverse each pixel point in the target region, calculate the difference between the pixel point coordinates and the belonging geometric center coordinates, and obtain the original displacement vector containing horizontal and vertical components. For the pixel points in the background region, set their vector values to zero.

[0111] Unitization processing: In order to eliminate the scale difference caused by different size target regions, the original displacement vector is normalized. The module length of the vector, that is, the Euclidean distance from the pixel to the center, is calculated, and the horizontal component and the vertical component are divided by the module length respectively to generate a standardized centripetal unit vector field label, which is used to calculate the training loss.

[0112] Regression prediction mechanism: At the end of the three-dimensional boundary reconstruction unit in the digital image processing architecture, a special pixel-level regression prediction module is designed:

[0113] Module structure: The module is composed of a point convolution layer, and the output channels correspond to the horizontal coordinate mapping and the vertical coordinate mapping of the vector field respectively; the network output layer does not use the activation function limiting the numerical range, so as to allow the output of negative values.

[0114] Numerical prediction: When the image to be processed is input, the module performs numerical regression on each spatial position of the feature map to directly predict the horizontal offset value and the vertical offset value of the point pointing to the center.

[0115] Field construction: The original offset value output by the network is converted into a centripetal unit vector field with clear direction and uniform module length through a pixel-by-pixel numerical normalization layer. The subsequent divergence calculation and Hough voting are based on the unit vector field.

[0116] It should be noted that when calculating the vector field regression loss, such as smooth L1 loss, only the pixel points (foreground pixels) inside the target region mask are calculated and back propagated; the pixel points in the background region are given zero weight or masked in the loss calculation to avoid gradient calculation errors caused by the inability to normalize zero vectors in the background.

[0117] The present application improves the robustness of adhesion segmentation by predicting the vector pointing to the geometric center through regression prediction instead of directly predicting the boundary, and has extremely high stability. Because in biological microscopic images, the edges of organoids are often blurred or overlapped, but the position of the geometric center (centroid) is usually determined and located inside the target. Let all pixels look at the centroid to form a centripetal force field with extremely high consistency, even if the edge is blurred, the attribution of the target can be accurately inferred through the vector flow, thereby effectively solving the dense adhesion problem.

[0118] Preferably, the texture gray scale analysis unit identifies the gray scale abnormal region in the following process:

[0119] The texture gray scale analysis unit receives the multi-scale filtering response map output by the image filtering unit; performs a multi-scale feature fusion operation to splice the features of the shallow layer containing high-resolution details and the features of the deep layer containing semantic information, and extracts features sensitive to local optical density changes through a spatial filtering matrix to obtain texture difference features representing the difference between the internal gray scale aggregation area and the surrounding gray scale smooth area of the target region;

[0120] Based on the extracted texture difference features, a pixel-level numerical mapping is performed to output a high-density region distribution map, wherein the pixel value represents the confidence of belonging to the gray abnormal region;

[0121] The target region mask is used as a spatial constraint condition to perform a spatial intersection operation with the high-density region distribution map to filter out noise pixels located outside the individual contour of the target region, ensuring that the high-density core is contained within the target region in terms of spatial topology;

[0122] The high-density region distribution map after topological filtering is subjected to Gaussian smoothing processing to fit the fuzzy boundary characteristics of the gray abnormal region, and an adaptive threshold is set for binary processing to generate a high-density core mask.

[0123] The spatial filtering matrix is specifically a learnable convolution kernel in a deep learning network, preferably a dilated convolution kernel with a dilated rate of 2 or 4, so as to automatically learn and extract feature response parameters for high light density texture regions during the training process.

[0124] The texture gray scale analysis unit performs multi-scale feature fusion to splice the high-resolution features of the shallow layer (retaining the jagged details of the necrotic region edge) and the semantic features of the deep layer (confirming that this is cell tissue rather than impurities). Through a spatial filtering matrix (convolution kernel), a response sensitive to local optical density changes is extracted. Generally, necrotic regions exhibit high optical density (black patches) and loss of texture uniformity.

[0125] First, a pixel-level high-density region distribution map is output, representing the probability of each pixel being necrotic; to remove interference such as stains and scratches in the background, a spatial topological constraint is introduced;

[0126] Operation logic: the generated target region mask is used as a spatial constraint condition to perform a spatial intersection operation with the high-density region distribution map. Based on biological facts, necrosis must occur within cell clusters. Therefore, any high-probability necrotic point located outside the target contour must be a false artifact and should be forcibly filtered out.

[0127] Post-processing: the filtered probability map is subjected to Gaussian smoothing processing to fit the fuzzy boundary characteristics of the natural diffusion of the necrotic region, and finally an adaptive threshold processing is performed to generate a binary high-density core mask.

[0128] The application solves the problem of poor noise resistance of the traditional threshold segmentation method by introducing a spatial topological constraint. In microscopic imaging, the scratch at the bottom of the culture dish and the edge of the bubble often have similar gray scale characteristics as the necrotic area. By using the determined target contour as prior knowledge, a logical constraint for spatial operation is constructed, greatly reducing the false positive detection rate and ensuring the biological reasonableness of the analysis results.

[0129] Further, the process of setting the adaptive threshold for the binarization processing specifically includes:

[0130] Extracting instance-level statistical features: traversing each independent target region mask, extracting the confidence values of all valid pixel points inside the target region mask in the high-density region distribution map, and calculating the confidence statistical mean and confidence standard deviation of the valid pixel points;

[0131] Constructing an instance-level dynamic threshold: based on the confidence statistical mean, a baseline is constructed, and the confidence standard deviation is multiplied by a preset sensitivity adjustment coefficient as a floating range. The baseline and the floating range are added to calculate the instance-level dynamic threshold for the target region;

[0132] Performing differential binarization: using the instance-level dynamic threshold, performing pixel-level comparison inside the corresponding target region in the high-density region distribution map, and marking the pixel points with confidence values higher than the instance-level dynamic threshold as high-density cores, so that different segmentation standards are applied to target regions with different contrasts in the same image.

[0133] The traditional global fixed threshold method is prone to serious over-segmentation (misjudging healthy tissue as necrotic) or under-segmentation (missing shallow necrosis) when facing batch-to-batch staining differences or microscope light source fluctuations.

[0134] The present scheme adopts a dynamic strategy based on instance context:

[0135] Instance-level statistical feature extraction: instead of macroscopic analysis of the entire image, individual organoids are analyzed one by one. Using the generated target region mask, the image is divided into several independent regions of interest. For the Nth organoid, only the high-density region distribution image pixel values within its mask range are extracted. Statistical analysis is performed on these pixel values to calculate two key indicators:

[0136] Confidence statistical mean: represents the average necrosis prediction probability inside the organoid, reflecting the overall suspected necrosis degree of the organoid.

[0137] Confidence standard deviation: represents the contrast difference of the texture inside the organoid. If the standard deviation is large, the core and edge boundary are distinct. If the standard deviation is small, the internal texture is uniform.

[0138] Construction of instance-level dynamic threshold: a unique threshold is calculated for each organoid; the calculation logic is: threshold = confidence statistical mean + (sensitivity adjustment coefficient x confidence standard deviation); wherein the sensitivity adjustment coefficient is a preset empirical constant, for example, the value is between 0.5 and 1; what is sought is a pixel point that is more significant than the average level and whose significant degree exceeds a certain fluctuation amplitude.

[0139] Differential binarization processing: using the calculated dynamic threshold, the inside of the specific organoid is binarized;

[0140] Scenario A, high contrast: for a necrotic core that is deeply stained and has clear boundaries, the standard deviation is large, and the calculated threshold will automatically increase, thereby accurately outlining the core edge and avoiding selecting the surrounding fuzzy transition zone.

[0141] Scenario B, low contrast / exposure insufficient: for an image that is overall dark and the network prediction confidence is generally low (for example, the maximum value is only 0.4), the mean value is also low. At this time, the calculated threshold will automatically decrease, for example, to 0.3, so that the relatively darkest core area can still be sensitively captured, avoiding missed detection.

[0142] Through the above steps, all organoids in the whole image are traversed, and a high-density core mask containing all individual internal details is generated.

[0143] The present application overcomes the interference caused by imaging quality fluctuations. In biological experiments, the brightness and contrast of microscopic images are greatly affected by staining time, reagent batch, and light source intensity. Traditional fixed threshold methods often can only adapt to images under specific conditions. The adaptive strategy proposed in this scheme takes the organoid itself as the reference system and calculates the relative statistical threshold. Whether the image as a whole is bright or dark, as long as there is a gray / texture difference between the necrotic core and healthy tissue, this algorithm can accurately separate them, greatly improving the robustness and environmental adaptability of the system, solving the problem of missed detection of weak feature samples, and significantly improving the detection rate of early necrotic cores.

[0144] Preferably, the process of constructing a three-dimensional geometric topology map comprises:

[0145] The generated target region mask is labeled with a connected domain, and each independent connected domain is instantiated as a grid node to construct a set of grid nodes.

[0146] Based on the target region mask, the geometric feature of each grid node is calculated, which at least includes the equivalent diameter and the roundness; meanwhile, based on the high-density core mask, the pixel statistical feature of each node is calculated, which at least includes the area proportion of abnormal gray region and the average gray of abnormal gray region; the geometric feature and the pixel statistical feature are combined to generate the initial feature vector of each node;

[0147] The Euclidean distance between the centroids of the target regions corresponding to any two grid nodes is calculated; if the Euclidean distance is less than the preset interaction threshold, it is determined that there is spatial interaction between the two, and a connection edge is established between the corresponding two nodes, thereby generating a three-dimensional geometric topology graph containing a node set, an edge set and a node feature vector.

[0148] Node instantiation and feature extraction: the target region mask is subjected to connected component analysis, and each independent organoid instance is instantiated as a grid node.

[0149] The attributes of each node are calculated:

[0150] Geometric feature: including equivalent diameter (area corresponding to the diameter of a circle), roundness (ratio of circumference to area).

[0151] Pixel statistical feature: based on the high-density core mask, the area proportion of necrotic region and the average gray value of necrotic region are calculated.

[0152] Organoids usually present spherical or ellipsoidal morphology in suspension culture or matrigel; although the image is a two-dimensional projection, there is a strong correlation between the projected area and the physical radius. Based on the projected area, the physical radius of the organoid is estimated, and it is assumed that the organoid is a sphere with this radius as a parameter in space.

[0153] Spatial interaction judgment distance calculation: the Euclidean distance between the centroids of any two nodes is calculated;

[0154] Interaction threshold: this threshold is not a fixed value, but is set based on the sum of the physical radii of the two, including the allowed small deformation tolerance; for organoids settled in the same focal plane, i.e. single-layer culture, if the centroid distance is less than the sum of the radii, it means that the two spheres have contacted or been squeezed in physical space; for three-dimensional culture, the out-of-focus targets are filtered out by preprocessing, and only the targets in the same focal plane are analyzed.

[0155] The connection edge is established between the nodes that meet the above conditions, thereby constructing a spatial geometric topology graph.

[0156] Further, the embodiment is preferably applicable to the scenario of suspension culture or matrix gel embedding culture, in which the organoids are affected by liquid surface tension or isotropic pressure of the matrix gel, and their morphology naturally tends to be a sphere or an ellipsoid, so the physical radius estimated based on the two-dimensional projection area has a high confidence.

[0157] The application breaks through the limitations of two-dimensional image analysis and restores the physical occupation relationship of organoids in the real culture environment through the "sphere assumption". Traditional two-dimensional distance calculation often ignores the influence of object size, and large objects are more likely to contact than small objects. The constructed topology graph truly reflects the relationship and provides a key structured data basis for subsequent differentiation of necrosis caused by physical crowding.

[0158] Preferably, the process of processing the three-dimensional geometric topology graph using a numerical iterative aggregation algorithm includes:

[0159] Input the three-dimensional geometric topology graph into the three-dimensional topology calculation module, calculate the data similarity between each target node and all connected neighbor nodes, generate normalized correlation coefficients to represent the influence weight of neighbor nodes on the three-dimensional space state of the target node, and generate the aggregation feature fused with the environmental context information.

[0160] Based on the correlation coefficient, the initial feature vector of the neighbor node is weighted and summed to transfer the morphological crowding information and gray distribution information of the neighbor to the target node, and generate the aggregation feature fused with the environmental context information.

[0161] Input the aggregation feature into the spatial mapping matrix for feature transformation, update the feature vector of the target node, and obtain the three-dimensional spatial morphological feature.

[0162] Preferably, the process of analyzing the three-dimensional spatial morphological feature includes:

[0163] Traverse all nodes in the three-dimensional geometric topology graph, calculate the local physical crowding degree of each target region based on the number of connection edges of each node combined with the distance weight sum of neighbor nodes.

[0164] The updated node feature vector, local physical crowding degree and gray abnormal area area ratio are used as joint input features; through feature space clustering analysis, the distribution difference of different gray modes in the feature space is mined to judge the state type of the target region; and the image feature parameters are obtained through the geometric morphological feature and state type of the grid node.

[0165] The architecture of the three-dimensional topology calculation module includes a geometric constraint calculation unit, a feature correlation calculation unit and a multi-modal fusion unit; the process of calculating the data similarity of all connected neighbor nodes and generating normalized correlation coefficients includes:

[0166] Computing geometric affinity weight: through the geometric constraint computing unit, the spatial Euclidean distance of the target node and the neighborhood node in the three-dimensional geometric topology graph is calculated, and the spatial Euclidean distance is mapped into a geometric affinity weight by using a Gaussian kernel function, so as to represent the adjacency strength in the physical space;

[0167] Computing semantic similarity weight: through the feature correlation computing unit, the dot product or cosine similarity between the initial feature vector of the target node and the initial feature vector of the neighborhood node is calculated, a semantic similarity weight is generated, and the correlation degree in the biological trait is represented;

[0168] Performing double-flow coupling normalization: through the multi-modal fusion unit, the geometric affinity weight and the semantic similarity weight are subjected to element-by-element multiplication or weighted splicing operation, a coupling weight is generated, and the coupling weight is subjected to Softmax normalization processing to obtain the final normalized correlation coefficient.

[0169] The traditional graph algorithm often only focuses on topological connection or only focuses on feature similarity. The scheme designs a "physical-pathological" double-factor coupling computing architecture.

[0170] Module architecture composition: the module is not a single black box model, but is composed of three parallel processing units:

[0171] Geometric constraint computing unit: responsible for processing the spatial position relationship between nodes, simulating the physical extrusion effect;

[0172] Feature correlation computing unit: responsible for processing the image attribute relationship between nodes, simulating the signal conduction effect;

[0173] Multi-modal fusion unit: responsible for combining physical information and visual information to generate the final attention distribution.

[0174] Calculation of geometric affinity weight: for each connection edge in the topology graph, i.e. connecting the target node and the neighborhood node, the centroid coordinates of the two are first extracted;

[0175] Distance measurement: the Euclidean distance of the two in the fitted three-dimensional space is calculated;

[0176] Gaussian mapping: since the closer the distance, the greater the influence, and the influence decays nonlinearly with the distance, the system uses a Gaussian kernel function to nonlinearly map the distance. The specific operation is as follows: the square of the distance is calculated, divided by a preset bandwidth parameter to control the decay rate, and then the negative number is taken to the exponential. The value obtained by this is the geometric affinity weight. The value ranges from 0 to 1, and the closer the distance, the closer the value to 1.

[0177] Semantic similarity weight calculation: Extract the initial feature vector of the target node and its neighborhood nodes, including diameter, roundness, necrosis rate, etc. First, project the feature vector to a common feature subspace through a linear mapping matrix to align the feature dimensions. Calculate the dot product (inner product) of the two projected vectors. The dot product result reflects the degree of similarity between the two organoids in morphology and pathological state, for example: both are large and necrotic organoids, the dot product value is large. This value is the semantic similarity weight.

[0178] Dual-flow coupling and normalization: In order to comprehensively consider the two factors of "close" and "look like", the multi-modal fusion unit performs a coupling operation: multiply the corresponding geometric closeness weight and semantic similarity weight.

[0179] This multiplication mechanism plays a "gating" role: only when two nodes are in close physical contact and highly related in state, their interaction will be significantly amplified.

[0180] Softmax normalization: For all neighbors of the target node, collect the coupling weights calculated above and perform Softmax normalization; that is, make the sum of all neighbor weights equal to 1, and the resulting value is the final normalized correlation coefficient.

[0181] Using the normalized correlation coefficient, the feature vectors of the neighbors are weighted and summed to complete a numerical iteration. This process can be repeated and stacked multiple times, allowing the node to perceive the state of "neighbors of neighbors" and obtain an updated target node feature vector, i.e., a three-dimensional spatial morphology feature.

[0182] The reason for using multiplication in dual-flow coupling is: compared to using addition, if two nodes are physically isolated but highly similar in biological state, the multiplication result tends to 0, effectively preventing information transfer between unrelated nodes; on the contrary, only when two nodes are physically adjacent and similar in state, the product is significant, this "and" logic is more consistent with the biological fact that crowding stress is locally concentrated. Softmax normalization converts the coupling weight into a probability distribution, making the aggregation operation on neighbors have a probabilistic explanation of attention, while ensuring that the weight sum is 1, preventing exponential growth of feature amplitude in multi-layer stacking.

[0183] The present application forces the model to pay attention to the physical distance by introducing the geometric affinity weight. In drug screening, the necrosis caused by environmental stress often has local aggregation, and the double-flow coupling mechanism can accurately capture this physical neighborhood-based collaborative change, thereby more accurately identifying the crowded stress pattern; by introducing the semantic similarity weight, it has adaptive selection capability, so that it can automatically ignore the noise nodes around when aggregating information, and only pay attention to the effective neighbors with reference value, thereby outputting purer and more discriminative image feature parameters.

[0184] The three-dimensional spatial morphology feature analysis is based on the updated feature vector and performs multi-dimensional feature analysis:

[0185] Local physical crowding degree calculation: based on the number of connected edges of the node and the neighborhood distance weight calculation.

[0186] Feature space clustering analysis: map the data to a high-dimensional feature space, and use clustering algorithms (such as K-Means or spectral clustering) to mine distribution patterns.

[0187] Pattern one, crowded stress: high physical crowding degree is strongly correlated with high necrosis rate, and the necrosis area tends to be distributed on the pressing contact surface.

[0188] Pattern two, drug toxicity: the correlation between necrosis rate and crowding degree is weak, and the necrosis phenomenon is uniformly distributed in the whole well plate and does not depend on the local density.

[0189] According to the clustering results, the state type of each organoid is determined, and the final image feature parameters are output; and the image is automatically labeled with information through the image feature parameters.

[0190] The present application uses a numerical iterative aggregation algorithm to realize the leap from individual visual features to group ecological features. By aggregating environmental context information, the algorithm can effectively distinguish between physical damage caused by pressing and chemical damage caused by drugs. This has very high value in drug screening and can help researchers exclude false positive results caused by too high culture density and accurately screen out truly effective drug compounds.

[0191] Example two:

[0192] The present application proposes an organoid automatic labeling method based on machine vision and deep learning, comprising:

[0193] Obtain the image to be processed under the high-throughput culture environment, and perform image preprocessing;

[0194] Input the preprocessed image to be processed into a preset digital image processing architecture; the digital image processing architecture includes an image filtering unit, a three-dimensional boundary reconstruction unit, and a texture gray scale analysis unit.

[0195] By the three-dimensional boundary reconstruction unit, pixel clustering based on a gradient vector field is performed, a centripetal unit vector field of each pixel point in the image pointing to a local geometric center to which the pixel point belongs is calculated, and dense adhesion between connected domains is numerically segmented based on a divergence characteristic of the vector field to generate an independent target region mask;

[0196] By the texture gray scale analysis unit, a pixel gray scale statistical analysis task of the image is performed; under a topological inclusion constraint, a gray scale abnormal region of the target region is extracted to generate a high-density core mask;

[0197] Based on the target region mask and the high-density core mask, a three-dimensional geometric topological graph of the organoid is constructed; a three-dimensional spatial morphological characteristic of a neighborhood node is obtained by using a numerical iterative aggregation algorithm; and the three-dimensional spatial morphological characteristic is analyzed to output an image feature parameter.

[0198] Specifically, the image acquisition and preprocessing process includes:

[0199] The input RGB color image (8-bit unsigned integer per channel, value range [0, 255]) is converted to the Lab color space, where the L channel represents the perceived brightness. The L channel is extracted as a single-channel gray scale image.

[0200] The CLAHE (Contrast Limited Adaptive Histogram Equalization) algorithm is applied to the L channel image, and the processing parameters are set as follows: the image segmentation block size is 8x8 pixels, a 4-pixel overlap is used between adjacent blocks to form a smooth transition band, and the histogram clipping limit value is set to 2.0. The specific process is as follows: the image is divided into non-overlapping 8x8 pixel basic blocks, the histogram of each block is calculated and equalized independently; at the block boundary, a bilinear interpolation algorithm is used for fusion, and the interpolation kernel size is 4x4 pixels to eliminate the seam artifacts between blocks.

[0201] The pixel values of the output image obtained after CLAHE processing are remapped to the [0, 1] floating point range to be sent to the subsequent neural network module for processing; the purpose of this preprocessing step is to suppress the effects of noise and uneven illumination commonly seen in microscopic imaging while maintaining the sharpness of the target edge gradient, thereby improving the accuracy of the subsequent segmentation module.

[0202] The design of the digital image processing architecture includes:

[0203] The image filtering unit fuses low-frequency global contour features and high-frequency local detail features, adapts to the scale changes of the target region at different stages, and outputs a multi-scale filtering response map;

[0204] The three-dimensional boundary reconstruction unit processes gradient edges and low-frequency contour information in the image to construct geometric shapes;

[0205] The texture gray scale analysis unit processes high-frequency texture details in the image and gray scale statistical characteristics, calculates and separates gray scale differences of different pixel density regions.

[0206] The image filtering unit adopts a deep residual network (ResNet-50) as a backbone network, network initialization adopts a Kaiming normal distribution initialization method, and weights are sampled from a normal distribution with a mean of 0 and a standard deviation of , where N is the number of input channels. If a pre-trained model is used, the weights pre-trained on the ImageNet dataset are used for parameter initialization, and fine-tuning is performed on the organoid dataset of the application.

[0207] The network input image resolution is uniform at 512x512 pixels, contains a single L channel (gray scale), the data type is float32, and the pixel value range is [0, 1].

[0208] The ResNet-50 backbone network includes 4 residual stages (Stage1-4), each stage generates feature maps C1, C2, C3, and C4 in turn, with resolutions of 1 / 4, 1 / 8, 1 / 16, and 1 / 32 of the original input, and channel numbers of 256, 512, 1024, and 2048, respectively. Each residual block contains 2 or 3 convolutional layers (depending on the block type), and each convolutional layer is configured with a batch normalization (BN) layer and a ReLU activation function; the residual connection adopts a bottleneck design with channel dimension reduction, with the number of channels in the middle layer being 1 / 4 of the input channel number of the block to reduce computational complexity.

[0209] A feature pyramid structure is constructed:

[0210] Features are extracted from C4, C3, C2, and C1, denoted as P5, P4, P3, and P2; the top-down path starts from P5 and performs 2x bilinear interpolation upsampling layer by layer to increase the spatial resolution to the size of the next layer; after each upsampling, the channel number is uniformly adjusted to 256 through a 1x1 convolution; the upsampled features are added element by element with the bottom features (C3, C2, C1) of the same level; after fusion, a 3x3 convolution layer is used for smoothing, and finally P5, P4, P3, and P2 four multi-scale feature maps are output, with channel numbers of 256.

[0211] Specifically, the image filtering unit outputs four feature maps of different scales, denoted as P2, P3, P4, and P5, with decreasing resolutions. In order to enable the three-dimensional boundary reconstruction unit and the texture gray scale analysis unit to simultaneously utilize global semantics and local details, this embodiment adopts a multi-level feature cascade fusion strategy;

[0212] Feature alignment: upsample the low-resolution feature maps P3, P4, P5 by bilinear interpolation algorithm to make their spatial dimensions consistent with the highest resolution feature map P2.

[0213] Channel concatenation: concatenate the feature maps of the four scales in the channel dimension. For example, if the channel number of P2 to P5 is 256, the total channel number of the concatenated feature map is 256x4, i.e. 1024.

[0214] Feature compression and redistribution: input the concatenated total feature map into a 3*3 convolution layer (including batch normalization BN and ReLU activation) to reduce the channel number to 256. The fused feature map is used as a common input and is sent to the decoding head of the three-dimensional boundary reconstruction unit and the texture gray analysis unit, respectively.

[0215] Further, the process of constructing a three-dimensional geometric topology graph includes:

[0216] Label the generated target region mask with connected domain, and instantiate each independent connected domain as a grid node to construct a set of grid nodes.

[0217] Based on the target region mask, calculate the geometric morphological features of each grid node, including at least the equivalent diameter and the circularity; at the same time, based on the high-density core mask, calculate the pixel statistical features of each node, including at least the area proportion of gray abnormal area and the average gray of gray abnormal area; combine the geometric morphological features and the pixel statistical features to generate an initial feature vector for each node.

[0218] Calculate the Euclidean distance between the centroids of any two grid nodes corresponding to the target regions; if the Euclidean distance is less than a preset interaction threshold, it is determined that there is a spatial interaction between the two, and a connection edge is established between the corresponding two nodes, thereby generating a three-dimensional geometric topology graph containing a node set, an edge set and a node feature vector.

[0219] The formula for determining the spatial relationship according to the interaction threshold is:

[0220] ;

[0221] wherein, is the Euclidean distance between the centroids i and j of the two target regions; , is the equivalent physical radius estimated based on the projected area; is a small deformation tolerance, to compensate for the non-rigid deformation of the organoid as a soft tissue when extruded; for example, the present embodiment takes Defined as 10% of the average of all equivalent physical radii; in addition, it can also be set to 3 to 10 pixel units based on the conventional microscopic imaging resolution.

[0222] when When the sum of the radii is less than the sum of the radii, it indicates that physical compression has occurred. Furthermore, when... When the sum of the radii is greater than the sum of the radii but less than the threshold, it indicates that the two are in a critical state of close adjacency but without significant deformation; both of these cases are considered to have spatial interaction and require the establishment of topological connection edges.

[0223] Preferably, the process of processing the 3D geometric topology map using a numerical iterative aggregation algorithm includes:

[0224] The three-dimensional geometric topology graph is input into the three-dimensional topology calculation module. For each target node, the data similarity between it and all connected neighboring nodes is calculated, and a normalized correlation coefficient is generated to characterize the influence weight of neighboring nodes on the three-dimensional spatial state of the target node.

[0225] Based on the correlation coefficient, a weighted summation operation is performed on the initial feature vectors of the neighboring nodes, and the morphological crowding information and gray-scale distribution information of the neighborhood are transmitted to the target node to generate aggregated features that integrate environmental context information.

[0226] The aggregated features are input into the spatial mapping matrix for feature transformation, and the feature vector of the target node is updated to obtain the three-dimensional spatial morphological features.

[0227] The spatial mapping matrix is ​​essentially a learnable linear transformation layer, containing a weight matrix and a bias vector. In this embodiment, the linear transformation layer preferably uses the LeakyReLU function with a slope of 0.2 on the negative half-axis.

[0228] Simple linear mappings cannot fit complex biological feature distributions. By introducing bias vectors and nonlinear activation functions, the model can perform nonlinear manifold mappings of crowded and pathological states in a high-dimensional feature space. This allows the model to distinguish between environmental stress and drug toxicity through simple cluster analysis in subsequent steps.

[0229] Preferably, the process of analyzing three-dimensional spatial morphological features includes:

[0230] Traverse all nodes in the 3D geometric topology graph, and calculate the local physical congestion of each target region based on the number of connecting edges of each node and the distance weights of neighboring nodes.

[0231] The updated node feature vector, local physical crowding degree and gray abnormal area area proportion are taken as joint input features; through feature space clustering analysis, the distribution difference of different gray modes in the feature space is mined to judge the state type of the target region; and the image feature parameters are obtained through the geometric morphological features of the grid nodes and the state type.

[0232] Specifically, the automatic state determination logic is: performing feature analysis on the multiple cluster centers generated by clustering. If the feature vector of a cluster center shows high local physical crowding degree, for example, the normalized value is greater than 0.7, and high necrosis rate, and the necrosis area tends to be distributed on the extrusion contact surface, then all the targets in the cluster are automatically marked as crowded stress state; if the feature vector of a cluster center shows low local physical crowding degree but high necrosis rate, then all the targets in the cluster are automatically marked as drug toxicity state; the rest of the clusters are marked as healthy state.

[0233] Although embodiments of the present application have been shown and described, it is to be understood that various modifications, substitutions, replacements and changes can be made to these embodiments without departing from the principles and spirit of the present application, and the scope of the present application is defined by the appended claims and their equivalents.

Claims

1. An automatic organoid annotation method based on machine vision and deep learning, characterized in that, include: Acquire images to be processed in a high-throughput culture environment and perform image preprocessing; The preprocessed image to be processed is input into a preset digital image processing architecture; the digital image processing architecture includes an image filtering unit, a three-dimensional boundary reconstruction unit, and a texture grayscale analysis unit. The image filtering unit integrates low-frequency global contour features and high-frequency local detail features to adapt to the scale changes of the target region at different stages and outputs a multi-scale filtering response map. The three-dimensional boundary reconstruction unit processes gradient edges and low-frequency contour information in the image to construct geometric shapes; The texture grayscale analysis unit processes high-frequency texture details and grayscale statistical features in the image, and calculates and separates the grayscale differences in regions with different pixel densities. Through the three-dimensional boundary reconstruction unit, pixel clustering based on gradient vector field is performed to calculate the centripetal unit vector field of each pixel in the image pointing to its local geometric center. Based on the divergence feature of the vector field, numerical segmentation is performed on the dense adhesion between connected domains to generate independent target region masks. The process of obtaining the target region mask includes: Calculate the normalized direction vector of each pixel in the image pointing to its local geometric center, forming a centripetal unit vector field; at the same time, calculate the membership degree of each pixel to the target foreground region, forming a foreground probability distribution map. Divergence operation is performed on the centripetal unit vector field to generate a vector field divergence map; the divergence operator is used to detect abrupt change regions in the vector direction field, and the positive divergence peak region is defined as the adhesion boundary barrier potential energy between target regions. Hough voting is performed on the centripetal unit vector field to count the convergence center of the vector flow. Regions with convergence intensity exceeding a preset threshold are identified as geometric center markers and used as water injection seeds for the watershed algorithm. Starting from the geometric center mark, a ground energy surface is constructed by superimposing the vector field divergence map onto the foreground probability distribution map, and a water-spreading operation is performed to generate mutually independent target region masks. The texture grayscale analysis unit performs pixel grayscale statistical analysis of the image; under the topological inclusion constraint, it extracts grayscale anomaly regions of the target region and generates a high-density core mask. Based on the target region mask and the high-density core mask, a three-dimensional geometric topology map of organoids is constructed; a numerical iterative aggregation algorithm is used to process the map to obtain the three-dimensional spatial morphological features of the neighboring nodes; based on the three-dimensional spatial morphological features, the image feature parameters are output.

2. The automatic organ-likeness annotation method based on machine vision and deep learning according to claim 1, characterized in that: The process of obtaining the geological energy surface includes: The foreground probability distribution map is numerically inverted to construct a basic topographic map, so that the interior area of ​​the target area is mapped as a low-potential basin, while the background area is mapped as a high-potential highland. Linear rectification is performed on the vector field divergence plot to filter out negative values ​​and retain positive value regions. The adhesion boundary features representing the vector divergence trend are extracted to construct the boundary barrier potential energy map. The base topographic map and the boundary barrier potential energy map are added pixel by pixel to generate the final geopotential surface.

3. The automatic organ-likeness annotation method based on machine vision and deep learning according to claim 1, characterized in that: The process by which the texture grayscale analysis unit identifies grayscale anomaly regions includes: The texture grayscale analysis unit receives the multi-scale filter response map output by the image filtering unit; performs a multi-scale feature fusion operation, splices the features containing high-resolution details in the shallow layer with the features containing semantic information in the deep layer, and extracts the feature response sensitive to local light density changes in the image through the spatial filtering matrix to obtain the texture difference features that characterize the difference between the grayscale cluster area inside the target area and the surrounding grayscale smooth area. Based on the extracted texture difference features, pixel-level numerical mapping is performed to output a high-density region distribution map, where the pixel value represents the confidence level of belonging to the gray-level abnormal region; The target region mask is used as a spatial constraint condition, and spatial intersection operation is performed with the high-density region distribution map to filter out noise pixels located outside the individual contours of the target region, ensuring that the high-density core is contained within the target region in spatial topology. Gaussian smoothing is applied to the high-density region distribution map after topological filtering to fit the fuzzy boundary characteristics of gray-level anomaly regions, and an adaptive threshold is set for binarization to generate a high-density core mask.

4. The automatic organ-likeness annotation method based on machine vision and deep learning according to claim 1, characterized in that: The process of constructing a 3D geometric topology graph includes: The generated target region mask is labeled with connected components, and each independent connected component is instantiated as a grid node to construct a set of grid nodes; Based on the target region mask, the geometric features of each grid node are calculated, including at least the equivalent diameter and roundness; simultaneously, based on the high-density core mask, the pixel statistical features of each node are calculated, including at least the area ratio of gray-level anomaly regions and the average gray level of gray-level anomaly regions; the geometric features and pixel statistical features are combined to generate the initial feature vector of each node. Calculate the Euclidean distance between the centroids of the target regions corresponding to any two grid nodes; if the Euclidean distance is less than a preset interaction threshold, it is determined that there is a spatial interaction between the two, and a connecting edge is established between the corresponding two nodes to generate a three-dimensional geometric topology graph containing a set of nodes, a set of edges, and node feature vectors.

5. The automatic organoid annotation method based on machine vision and deep learning according to claim 4, characterized in that: The process of processing 3D geometric topology graphs using numerical iterative aggregation algorithms includes: The three-dimensional geometric topology graph is input into the three-dimensional topology calculation module. For each target node, the data similarity between it and all connected neighboring nodes is calculated, and a normalized correlation coefficient is generated to characterize the influence weight of neighboring nodes on the three-dimensional spatial state of the target node. Based on the correlation coefficient, a weighted summation operation is performed on the initial feature vectors of the neighboring nodes, and the morphological crowding information and gray-scale distribution information of the neighborhood are transmitted to the target node to generate aggregated features that integrate environmental context information. The aggregated features are input into the spatial mapping matrix for feature transformation, and the feature vector of the target node is updated to obtain the three-dimensional spatial morphological features.

6. The automatic organ-likeness annotation method based on machine vision and deep learning according to claim 5, characterized in that: The process of analyzing three-dimensional spatial morphological features includes: Traverse all nodes in the 3D geometric topology graph, and calculate the local physical congestion of each target region based on the number of connecting edges of each node and the distance weights of neighboring nodes. The updated node feature vector, local physical congestion, and area ratio of gray-level anomalies are used as joint input features. Through feature space clustering analysis, the distribution differences of different gray-level modes in the feature space are explored to determine the state type of the target area. Image feature parameters are obtained through the geometric morphology features and state types of grid nodes.

Images