Organ-like drug sensitivity automatic detection method and system for multi-sample parallel processing

Abstract

The invention provides an organoid drug sensitivity automatic detection method and system for multi-sample parallel processing, and the method comprises the steps: obtaining the optical path length data of each hole site of a cell culture perforated plate, and determining the initial focal plane position of each hole site; matching the deformed feature point coordinates with the initial three-dimensional structure data of the similar organs, and calculating the focal plane offset caused by the deformation of the similar organs through least square fitting; adjusting the focal length of the optical imaging system according to the corrected focal plane position by adopting a fast focusing algorithm, and re-collecting an organ-like image; using a layer-by-layer stacked three-dimensional reconstruction algorithm to perform registration splicing on the recollected image sequence to obtain complete three-dimensional structure data after the organ-like body is deformed; and according to the complete three-dimensional structure data, extracting the volume, surface area and cell distribution density characteristics of the organoid, and generating a morphological characteristic report of the drug sensitivity experiment.

Description

Technical Field

[0001] The present invention relates to the field of information technology, and in particular, to an automatic organoid drug sensitivity detection method and system for parallel processing of multiple samples. Background Art

[0002] In the automatic organoid drug sensitivity detection for parallel processing of multiple samples, the performance of the optical imaging system directly determines the accuracy of organoid morphological feature extraction. The design of the multi-well plate distributes samples at different spatial positions, resulting in a contradiction between imaging clarity and depth of field range when the optical system acquires images. Specifically, there are slight differences in the optical path lengths of each well position of the multi-well plate, which may lead to insufficient image clarity at some well positions, while at other well positions, the three-dimensional structure of the organoid cannot be completely captured due to depth of field limitations. When the organoid deforms during the drug sensitivity experiment, the spatial distribution of its three-dimensional structure changes, such as changes in height, volume, or density, which further exacerbates the challenges of optical imaging. The focal plane tracking strategy is the key to solving this problem, but its adjustment method needs to match the characteristics of organoid deformation. Organoid deformation may cause key feature points on its surface or inside to deviate from the original focal plane. If the focal plane tracking strategy fails to respond in a timely manner, it may lead to blurred or distorted images in some key regions. This kind of blur or distortion will affect the accuracy of morphological feature extraction, such as the calculation of the volume, surface area, or cell distribution density of the organoid. Further, the adjustment of the focal plane tracking strategy needs to be combined with the optical characteristics of the multi-well plate. The material, pore size, and light transmittance of the multi-well plate will affect the refraction and scattering of light, thereby changing the actual position of the focal plane. Therefore, the focal plane tracking strategy not only needs to respond to the deformation of the organoid, but also needs to compensate for the optical interference of the multi-well plate itself. Summary of the Invention

[0003] The present invention provides an automatic organoid drug sensitivity detection method for parallel processing of multiple samples, mainly including:

[0004] Obtain the optical path length data of each well position of the cell culture multi-well plate, and determine the initial focal plane position of each well position;

[0005] According to the initial focal plane position, collect the organoid images of each well position, extract the organoid contour using an edge detection algorithm, and map the two-dimensional contour to a three-dimensional space model using a shape matching method to obtain the initial three-dimensional structure data of the organoid;

[0006] During the drug sensitivity experiment, according to the characteristics of the volume and surface shape changes of the organoid, use the optical flow method to track the position changes of the key feature points on the organoid surface to obtain the coordinates of the deformed feature points;

[0007] Match the coordinates of the deformed feature points with the initial three-dimensional structure data of the organoid, and calculate the focal plane offset caused by the deformation of the organoid through least squares fitting;

[0008] According to the focal plane offset, combined with the refractive index and absorption coefficient of the multi-well plate, use the ray tracing algorithm to simulate the refraction and scattering paths of light in the multi-well plate, and obtain the corrected focal plane position;

[0009] Adopt a fast focusing algorithm, adjust the focal length of the optical imaging system according to the corrected focal plane position, and re-acquire the organoid images;

[0010] Use a three-dimensional reconstruction algorithm of layer-by-layer stacking to register and splice the re-acquired image sequence, and obtain the complete three-dimensional structure data of the deformed organoid;

[0011] According to the complete three-dimensional structure data, extract the volume, surface area and cell distribution density characteristics of the organoid, and generate a morphological feature report for the drug sensitivity experiment.

[0012] Furthermore, obtain the optical path length data of each well position of the cell culture multi-well plate, and determine the initial focal plane position of each well position, including: collecting the optical path length data of each well position of the multi-well plate through a high-precision optical sensor, screening the collected data according to the preset upper and lower limits of light intensity values, calculating the optimal collection position of the optical sensor in the direction perpendicular to the plane of the multi-well plate using a Gaussian distribution, and obtaining the initial well position optical data set. Use an optical path sensor to scan the cell distribution in the multi-well plate layer by layer at a preset scanning interval, apply wavelet transform denoising to the scanned data, obtain a clarity evaluation index by calculating the variance of the gray value of each scanned image layer, and select the scanned layer corresponding to the maximum clarity index value for each well position as the focal plane candidate layer. Compare and calibrate the standard aperture size of the multi-well plate with the measured aperture size in the actual collected optical image, calculate the aperture ratio coefficient using the least squares method, and correct the spatial position of the focal plane candidate layer through this coefficient to obtain the well position focal plane spatial calibration data. For the optical data of each well position, establish a functional relationship between light intensity and focal plane offset based on a preset optical model, calculate the optical path stability coefficient using the maximum likelihood estimation method, and if the optical path stability coefficient is less than the preset threshold of 0.8, correct the focal plane position using the light intensity attenuation compensation function. Optimize the corrected focal plane position data through recursive least squares method, establish a mapping relationship between the optical path length and the focal plane position, and obtain the final initial focal plane position data of each well position.

[0013] Furthermore, according to the initial focal plane position, organoid images at each well position are collected, and the organoid contours are extracted using an edge detection algorithm. The two-dimensional contours are mapped to a three-dimensional spatial model using a shape matching method to obtain the initial three-dimensional structure data of the organoids, including: starting a layer-by-layer scan of the multi-well plate according to the initial focal plane position data, collecting a multi-level image sequence of the organoids at 5-micron intervals in the vertical direction for each well position, removing the acquisition noise through Gaussian mean filtering with a radius of 3 pixels, enhancing the image edge features using a 5×5 Laplacian operator to obtain an enhanced organoid image sequence. For the enhanced organoid image sequence, the organoid edge contour point sets are extracted using Sobel operators of three scales, 3×3, 5×5, and 7×7, the extracted contour point sets are subjected to non-maximum suppression processing, and a double-threshold method with a high threshold of 0.8 and a low threshold of 0.3 is used to connect the strong and weak edges of the contour point sets to obtain the edge contour feature data of the organoids. Polar coordinate transformation is used to perform spatial mapping on the contour feature data, a depth estimation function is established by calculating the gray gradient between adjacent layer images, the surface depth information of the organoids is calculated according to the correspondence between the gray gradient value and the depth, and the Delaunay triangulation algorithm is used to perform grid processing on the contour feature points to obtain the initial spatial grid structure data of the organoids. Curvature, volume, and surface area characteristic parameters of the organoids are extracted from the initial spatial grid structure data, the spatial reconstruction of the organoid morphological features is performed through a convolutional neural network, and the Laplacian coordinate-based grid smoothing algorithm is applied to the reconstruction result to optimize the surface morphology. An iterative grid deformation algorithm is used to perform local adjustment on the three-dimensional structure, the grid vertex displacement amount is calculated by the least squares method, and the grid morphology is corrected based on a preset deformation threshold to obtain the final three-dimensional structure data of the organoids.

[0014] Furthermore, during the drug sensitivity experiment, according to the volume and surface shape change characteristics of the organoids, the optical flow method is used to track the position changes of the key feature points on the surface of the organoids to obtain the coordinates of the deformed feature points, including: extracting the area with a curvature value greater than the preset threshold on the surface of the organoids as the feature point sampling position through the region growing method, calculating the local curvature value and gradient direction for the sampling position, screening the feature points using the Harris corner detector with a response value greater than the preset threshold to obtain the initial surface feature point data. Calculate the local gray gradient change direction and change amplitude based on the initial surface feature points, perform a four-layer pyramid decomposition on the collected image sequence, use the dense optical flow algorithm with a window size of 16 pixels to calculate the feature point motion vector field, and obtain the surface deformation feature data. Establish a Markov state transition matrix for the surface deformation feature data, predict the morphological change trend of the organoids through the Bayesian prediction algorithm, set the state prediction error threshold to screen the prediction results, and obtain the feature point position prediction data. Use a particle filter with 200 particles to correct the trajectory of the feature point position prediction data, perform state estimation by calculating the particle weights, set the minimum weight threshold to eliminate abnormal trajectories, and obtain the feature point motion trajectory data. For the feature point motion trajectory data, calculate the spatial coordinates of the feature points using three-dimensional orthogonal basis tensor projection, calculate the displacement of the surface grid nodes of the organoids through triangulation, and obtain the spatial coordinate data of the deformed feature points. Calculate the surface grid deformation amount of the organoids based on the spatial coordinate data of the deformed feature points, optimize the feature point trajectory through a regression neural network, set the maximum deformation threshold to correct abnormal deformations, and obtain the coordinates of the surface feature points of the deformed organoids after correction.

[0015] Further, match the coordinates of the deformed feature points with the initial three-dimensional structure data of the organoid, and calculate the focal plane offset caused by the deformation of the organoid by least squares fitting, including: using the nearest neighbor matching algorithm to pair the coordinates of the deformed feature points with the feature points in the initial three-dimensional structure data, establishing a point position correspondence relationship by calculating the Euclidean distance between the feature points, eliminating the point positions with a matching distance greater than three times the average distance between points to obtain the feature point matching data. Calculate the displacement of each point position in the three coordinate axis directions according to the feature point matching data, calculate the point position deformation weight by the weighted average method, screen the point positions with a weight value less than the preset threshold to obtain the effective deformation data. Construct a three-dimensional coordinate transformation matrix for the effective deformation data, extract the main transformation direction and transformation amplitude by singular value decomposition, retain the main transformation components by setting the singular value screening threshold to obtain the spatial transformation feature data. Based on the spatial transformation feature data, construct a five-layer deep neural network, where the number of nodes in the input layer is the same as the feature dimension, the number of nodes in the hidden layer is halved in turn, and the output layer corresponds to the displacement components in three directions of the focal plane to obtain the focal plane displacement prediction data. Apply principal component analysis to the focal plane displacement prediction data, select the eigenvectors with a cumulative contribution rate reaching the preset threshold as the main offset directions to obtain the focal plane offset feature data. Establish a least squares fitting equation set according to the focal plane offset feature data, use the conjugate gradient method to iteratively solve the equation set, and set the convergence threshold to control the iteration termination condition to obtain the focal plane offset data caused by the deformation of the organoid.

[0016] Furthermore, based on the focal plane offset, combined with the refractive index and absorption coefficient of the porous plate, the ray tracing algorithm is used to simulate the refraction and scattering paths of light in the porous plate, and the corrected focal plane position is obtained, including: constructing a three-dimensional grid structure of the porous plate according to the focal plane offset, collecting refractive index and absorption coefficient data for each grid cell, calculating the material parameter distribution between grid nodes by cubic spline interpolation method, evaluating the interpolation accuracy by cross-validation method, and obtaining the optical property distribution data of the porous plate. Constructing a Monte Carlo ray tracing model for the optical property distribution data of the porous plate, setting the number of rays as the number of incident rays at each hole position, generating the spatial distribution of incident rays by stratified sampling method, calculating the propagation direction of each ray in the porous plate by backward tracing, and obtaining the ray spatial propagation data. According to the ray spatial propagation data, Snell's law is used to calculate the refraction angle at each interface, the interface reflection loss is calculated by the Franck-Condon principle, and the interfaces with reflectivity greater than the preset threshold are marked to obtain the interface optical property data. Calculating the scattered light intensity distribution for the interface optical property data, using a deep neural network to spatially model the scattering loss, constructing a scattered intensity attenuation function by recursive calculation, and obtaining the light intensity attenuation distribution data. Establishing the Lambert-Beer attenuation equation according to the light intensity attenuation distribution data, determining the light intensity compensation coefficient at each spatial position by iterative calculation, and performing spatial smoothing on the compensation coefficient to obtain the light intensity compensation data. Constructing an optimization equation set for the focal plane position for the light intensity compensation data, solving the equation set by conjugate gradient method, setting a convergence threshold to control the number of iterations, and eliminating abnormal solutions by standard deviation analysis to obtain the corrected focal plane position data.

[0017] Furthermore, obtain the refractive index distribution data of the porous plate, calculate the wavefront distortion amount of the transmitted light, adjust the scattering parameters of the preset path random distribution model, fuse the path random distribution model with the absorption coefficient, derive the phase compensation data, and update the correction parameters of the focal plane offset amount, including: divide the side length into square grid units with a preset size according to the porous plate structure, measure the refractive index distribution for the grid units using a double-beam interferometer, extract the frequency components of the optical path difference through fast Fourier transform, set a cut-off frequency to filter out high-frequency noise, and obtain the refractive index distribution map of the porous plate. Perform wavefront distortion decomposition on the refractive index distribution map using the third-order Zernike polynomial, calculate the coefficients of each polynomial by the least squares method, set a coefficient threshold to screen the main distortion terms, and obtain the wavefront distortion eigenvector. Establish a three-layer neural network for the wavefront distortion eigenvector, with the input layer corresponding to the distortion coefficients, the hidden layer using the hyperbolic tangent activation function, and the output layer generating the scattering parameter distribution to obtain the updated path random distribution parameters. Establish a zero-order Bessel function based on the path random distribution parameters, calculate the coefficients of the higher-order terms through recurrence relations, perform spatial decomposition on the scattering intensity, and obtain the scattering intensity distribution function. Perform a convolution operation on the spatial distributions of the scattering intensity distribution function and the absorption coefficient, calculate the light intensity attenuation sequence through discrete sampling, generate a phase compensation matrix using singular value decomposition, and obtain the phase compensation data. Construct a focal plane offset mapping function based on the phase compensation data, use the conjugate gradient method to iteratively optimize the correction parameters, and set the gradient threshold and the maximum number of iterations as the convergence conditions to obtain the correction parameters of the focal plane offset amount.

[0018] Furthermore, a fast focusing algorithm is adopted to adjust the focal length of the optical imaging system according to the corrected focal plane position, and the organoid image is recollected, including: constructing a Gaussian beam propagation equation based on radial symmetry according to the corrected focal plane position, establishing a focal length adjustment curve by calculating the mapping relationship between the beam cross-sectional radius and the optical path length, optimizing the focal length adjustment curve by using the gradient descent method with adaptive learning rate to obtain the initial focal length adjustment parameters. An optical path compensation function based on the refractive index distribution is established for the initial focal length adjustment parameters, the optical path difference value is measured by a Michelson interferometer, the optical path difference data is decomposed by discrete wavelet transform at multiple scales, the reconstruction coefficient of the maximum energy scale is selected to obtain the optical path compensation data. A five-layer convolutional neural network is constructed according to the optical path compensation data, the input layer receives the optical path compensation vector, the middle three layers adopt the residual connection structure, the output layer generates the focal parameter curve, and the standardized focal curve is obtained through batch normalization. A focal length adjustment interval is generated for the standardized focal curve by using the uniform division method, the adjustment cost between adjacent intervals is described by the state transition matrix, and the minimum adjustment path is calculated by the forward dynamic programming algorithm to obtain the focusing adjustment sequence. Based on the focusing adjustment sequence, the piezoelectric ceramic actuator is controlled to adjust the lens group spacing, the image variance and gradient features are calculated by the Gaussian kernel support vector machine, the clarity of the collected image is evaluated to obtain the clarity evaluation data. A focal length fine-tuning function is established for the clarity evaluation data, the focal length deviation is calculated by the weighted least squares method, the focal length is adjusted in a closed loop by the proportional integral controller, and the dead zone threshold is set to control the adjustment accuracy to obtain the focal length parameters.

[0019] Furthermore, the layer-by-layer stacked three-dimensional reconstruction algorithm is used to register and splice the recollected image sequence to obtain the complete three-dimensional structure data of the deformed organoid, including: constructing a corner response function according to the recollected image sequence, extracting the feature point coordinates by setting the corner response threshold, screening the feature points by non-maximum suppression with a radius of the preset pixel value, calculating the local binary pattern feature by using a circular sampling window to obtain the image feature point data. A cross-correlation response function is constructed for the image feature point data, the image translation amount and rotation angle are calculated by the phase correlation method, the spectrum data is transformed by the fast Fourier transform, and the registration result is evaluated by combining the normalized cross-correlation coefficient to obtain the image registration data. The corresponding relationship of the feature points is calculated according to the image registration data, the corresponding points are screened by the feature point matching distance threshold, the inter-layer transformation matrix is calculated by the least squares method, and the singular value of the transformation matrix is decomposed to obtain the spatial transformation parameters. A five-layer deep neural network is constructed for the spatial transformation parameters, the image features are extracted by the convolutional layer, a four-layer spatial pyramid is constructed by setting the pooling kernel size, and the multi-scale features are fused by the skip connection to obtain the feature fusion data. The deconvolution reconstruction is carried out according to the feature fusion data, the image resolution is restored by bilinear interpolation, the structure contour is extracted by the edge detection operator, and the contour points are spatially registered to obtain the three-dimensional contour data.

[0020] Furthermore, based on the complete three-dimensional structure data, the volume, surface area, and cell distribution density characteristics of the organoids are extracted to generate a morphological feature report for the drug sensitivity experiment, including: According to the vertex coordinates and patch information in the complete three-dimensional structure data, a tetrahedral meshing function is constructed by the least squares criterion, the meshing size threshold is set to control the size of the tetrahedrons, and the three-point Gaussian integration formula is used to calculate the volume of the voxel units to obtain the initial volume data of the organoids. For the three-dimensional structure grid, the barycentric coordinate mapping is used to calculate the area of the triangular patches, a surface weighting function is constructed by the mean curvature and Gaussian curvature, and the curvature threshold is set to identify the surface wrinkled areas to obtain the initial surface area data. A five-layer convolutional neural network is constructed based on the initial surface area data, the feature data is processed by the batch normalization layer, the pooling layer is set to compress the feature dimension, and the fully connected layer is used to fuse the multi-scale features to obtain the surface area feature data. According to the three-dimensional structure gray scale map, the adaptive segmentation threshold is calculated by the maximum inter-class variance method, the three-dimensional connected regions are labeled by the twenty-six neighborhood criterion, the centroid calculation and statistics are performed on the connected regions to obtain the nucleus position data. A spatial density function is constructed for the nucleus position data, the local area density distribution is calculated by the kernel density estimation method, and the density gradient field is calculated by the radial basis support vector regression to obtain the cell density feature data. The volume data, surface area features, and cell density features are standardized, the top five morphological parameters with the highest feature contribution rates are extracted by principal component analysis, and a feature description function is constructed by cubic polynomial fitting to obtain the morphological feature report data.

[0021] The present invention provides an organoid drug sensitivity automatic detection system for parallel processing of multiple samples, mainly including:

[0022] An optical path length calculation module, which is used to obtain the optical path length data of each well position of the cell culture multi-well plate and calculate the initial focal plane position of each well position through a preset optical model;

[0023] An initial three-dimensional structure reconstruction module, which is used to collect the organoid images of each well position according to the initial focal plane position, extract the organoid contours by using the edge detection algorithm, and map the two-dimensional contours to the three-dimensional space model by using the shape matching method to obtain the initial three-dimensional structure data of the organoids;

[0024] A feature point tracking module, which is used to track the position changes of the key feature points on the surface of the organoids by using the optical flow method during the drug sensitivity experiment according to the volume and surface shape change characteristics of the organoids to obtain the coordinates of the deformed feature points;

[0025] A focal plane offset calculation module, which is used to match the coordinates of the deformed feature points with the initial three-dimensional structure data of the organoids and calculate the focal plane offset caused by the deformation of the organoids by least squares fitting;

[0026] A ray tracing simulation module, configured to simulate the refraction and scattering paths of light rays in a porous plate according to the focal plane offset, in combination with the refractive index and absorption coefficient of the porous plate, using the ray tracing algorithm to obtain the corrected focal plane position;

[0027] A focal length adjustment module, configured to use a fast focusing algorithm to adjust the focal length of the optical imaging system according to the corrected focal plane position and re-acquire the organoid image;

[0028] A three-dimensional reconstruction module, configured to use a layer-by-layer stacked three-dimensional reconstruction algorithm to register and splice the re-acquired image sequence to obtain the complete three-dimensional structure data of the deformed organoid;

[0029] A morphological feature generation module, configured to extract the volume, surface area, and cell distribution density features of the organoid according to the complete three-dimensional structure data and generate a morphological feature report for the drug sensitivity experiment.

[0030] The technical solution provided by the embodiment of the present invention may include the following beneficial effects:

[0031] The present invention discloses an automatic organoid drug sensitivity detection method for parallel processing of multiple samples. The method first calculates the initial focal plane position of each well in the porous plate through an optical model, acquires the organoid image and extracts its initial three-dimensional structure. During the experiment, the optical flow method is used to track the position changes of the surface feature points of the organoid, and the focal plane offset is calculated in combination with the initial structure data. Subsequently, the present invention considers the optical characteristics of the porous plate, uses the ray tracing algorithm to correct the focal plane position, and uses the fast focusing algorithm to adjust the focal length of the imaging system to re-acquire the image. Finally, the complete structure data of the deformed organoid is generated through the three-dimensional reconstruction algorithm, and the relevant features are extracted to generate a drug sensitivity experiment report. The present invention effectively solves the problem of imaging defocus caused by the deformation of the organoid in the drug sensitivity experiment, realizes the dynamic and accurate monitoring of the three-dimensional structure of the organoid, and provides reliable technical support for drug screening and effect evaluation. BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Figure 1 It is a flowchart of an automatic organoid drug sensitivity detection method for parallel processing of multiple samples according to the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0033] To further understand the content of the present invention, the present invention will be described in detail in combination with the drawings and embodiments. The following further describes the present application in detail with reference to the drawings and embodiments. It can be understood that the specific embodiments described herein are only used to explain the related invention, rather than limiting the invention. In addition, it should be noted that for the sake of description, only the parts related to the invention are shown in the drawings.

[0034] Such as Figure 1In this embodiment, a method for automatically detecting drug sensitivity of organoids by processing multiple samples in parallel may specifically include:

[0035] S101, collecting optical path length data of each well of a cell culture multi-well plate and calculating an initial focal plane position of each well using a preset optical model.

[0036] S1011. In practical applications, a high-precision optical sensor is used to measure the optical path length data of each well of a multi-well plate. Valid data is screened according to a preset light intensity range and an initial well optical data set is generated. A Gaussian distribution model is used to determine the optimal acquisition position of the sensor in a direction perpendicular to the multi-well plate to improve measurement stability. For example, a stable light intensity signal is obtained 15 mm above the culture well. Subsequently, the distribution of organoids in the multi-well plate is scanned layer by layer at 50-micron intervals by an optical path sensor. Wavelet transform is applied to the scanned data for denoising and the grayscale value variance of each layer of the image is calculated as a clarity index. The layer corresponding to the maximum clarity index is selected as the focal plane candidate layer to ensure that the candidate layer can reflect the optimal imaging plane of the organoid.

[0037] S1012. After determining the candidate focal plane layer, compare the standard aperture size of the porous plate with the aperture size measured in the actual optical image, calculate the aperture ratio coefficient by the least squares method, and use the coefficient to correct the spatial position of the candidate layer to obtain the focal plane spatial calibration data of the hole position; further establish a functional relationship between light intensity and focal plane offset for the calibration data, if the optical path stability coefficient calculated by the maximum likelihood estimation method is less than the preset threshold value of 0.8, use the light intensity attenuation compensation function to correct the focal plane position, and optimize the correction result by the recursive least squares method, establish a mapping relationship between the optical path length and the focal plane position, and generate the initial focal plane position data for each hole position with an accuracy of up to sub-micron level.

[0038] In the process of collecting optical path length data, the application of high-precision optical sensors is crucial. Taking a 96-well plate as an example, each well has a diameter of 6.4 mm and a depth of 11.2 mm. The optical sensor needs to ensure measurement accuracy to adapt to the slight differences between wells. Experiments have shown that the light intensity range is set to an upper limit of 850 units and a lower limit of 150 units, which can effectively eliminate abnormal data. Wavelet transform denoising uses the db4 basis function and performs a 4-layer decomposition to remove Gaussian white noise and ensure the reliability of the scanned data.

[0039] The selection of focal plane candidate layers is based on the calculation of grayscale value variance. Experiments show that the variance of the five layers of data near the focal plane is significantly higher than that of other areas. The layer corresponding to the maximum variance value is the best candidate layer. The calculation of the aperture ratio coefficient reflects the magnification characteristics of the optical system, which is usually between 0.95 and 1.05. If it deviates from 1.00 by more than 0.03, the equipment needs to be recalibrated.

[0040] The functional relationship between the light intensity and the focal plane offset shows an exponential decay law, which is affected by the refractive index of the culture medium and the cell density. If the optical path stability is insufficient, the compensation function uses the piecewise linear interpolation method, and different coefficients are applied for adjustment in different offset intervals. When optimizing by the recursive least squares method, the sliding window size is 10 data points, and the optimization converges when the mean square error is less than 0.1 micrometer. This method can effectively cope with the focal plane drift caused by temperature changes or mechanical vibrations, and is especially suitable for long-term drug sensitivity test scenarios.

[0041] It can be understood that the embodiments of the present invention do not overly limit the specific model or scanning interval of the optical sensor, which can be adjusted by technicians according to experimental requirements to ensure the accuracy and applicability of the focal plane position data. The obtained initial focal plane position lays a solid foundation for subsequent deformation tracking and three-dimensional reconstruction.

[0042] S102. Start the optical imaging system according to the initial focal plane position to collect the organoid images of each well position, extract the organoid contours using the edge detection algorithm, and construct the initial three-dimensional structure data through the shape matching method.

[0043] Based on the initial focal plane position, use the optical imaging system to scan each well position of the multi-well plate layer by layer to collect the multi-level image sequence of the organoid. The scanning interval is set to 5 micrometers to adapt to the subtle structural changes of the organoid with a diameter in the range of 200 to 800 micrometers. During the scanning process, the refractive index difference of the culture medium and the light intensity attenuation may introduce noise. Therefore, the image sequence is preprocessed by the Gaussian mean filter, and the filtering radius is set to 3 pixels, which can effectively smooth the random noise while retaining the key edge information of the organoid. Subsequently, the Laplace operator is applied to enhance the edge features of the image. The operator size is 5×5, the central weight is 4, and the surrounding weights are -1. This configuration can highlight the gray jump region and make the organoid contour clearer.

[0044] S1021. In the edge detection stage, use the multi-scale Sobel operator to extract the contour point set of the organoid. The operators with sizes of 3×3, 5×5, and 7×7 are respectively used to adapt to the edge features of different-shaped organoids. Among them, the 3×3 operator is suitable for detecting sharp edges, and the 7×7 operator is more suitable for extracting smooth edges; apply non-maximum suppression to the extracted contour point set to eliminate redundant edge points, and process the strong and weak edges through the double-threshold connection method. The high threshold is set to 0.8, and the low threshold is set to 0.3 to ensure edge continuity while reducing false edge interference, and generate the edge contour feature data of the organoid.

[0045] S1022. Based on the contour feature data, the two-dimensional edge point set is mapped into the three-dimensional space by polar coordinate transformation. A depth estimation function is established by analyzing the gray gradient between adjacent layer images. When the gradient value is greater than 50, it indicates a significant depth change. This function is used to calculate the depth distribution on the surface of the organoid. Subsequently, the contour feature points are meshed by the Delaunay triangulation algorithm. This algorithm generates high-quality triangular meshes based on the principle of maximizing the minimum angle. When the interior angle of a triangle is greater than 30 degrees, the grid stability can be ensured, and the initial spatial grid structure data is obtained, providing a basis for three-dimensional reconstruction.

[0046] For the initial spatial grid structure data, the morphological feature parameters of the organoid are further extracted. For example, the local curvature reflects the degree of surface wrinkles. Usually, a value between 0.02 and 0.05 indicates a normal development state. The volume is calculated by the grid closed surface integral. The daily growth rate of healthy organoids is generally between 15% and 25%. To optimize the grid quality, a grid smoothing algorithm based on Laplace coordinates is applied, with the weight coefficient set to 0.5 and the number of iterations controlled within 10 times, which can smooth the surface noise and retain the detailed features. In addition, an iterative grid deformation algorithm is used to locally adjust the three-dimensional structure. The displacement of the grid vertices is calculated by the least squares method, and the deformation threshold is set to 10% of the original grid side length to ensure the stability of the adjusted grid topology structure.

[0047] To improve the accuracy of three-dimensional reconstruction, a convolutional neural network can be introduced to perform spatial reconstruction on the morphological features of the organoid. The network is configured with 32 feature channels, and each channel extracts morphological information at different scales, such as edge sharpness, surface curvature, etc. High-dimensional feature maps are generated through multiple convolutional and pooling operations and then mapped back to the three-dimensional space. The reconstruction result is further optimized by combining with the Laplace smoothing algorithm to avoid surface irregular fluctuations caused by noise. The advantage of Delaunay triangulation is that the generated grid can adapt to the irregular shape of the organoid. Compared with the uniform grid, its computational efficiency is higher and it is more suitable for subsequent morphological analysis.

[0048] In practical applications, the double-threshold method of edge detection can effectively deal with the edge ambiguity of the organoid through the dynamic combination of high and low thresholds. For example, when the organoid deforms under the action of drugs, the edge may show a gradual change characteristic. An overly high single threshold may truncate the continuous edge, while an overly low threshold may introduce too much noise. The double-threshold method just balances the contradiction between the two. The gray gradient depth estimation relies on the light intensity change law between adjacent layers. Experiments show that when the refractive index of the culture medium is between 1.33 and 1.37, the linear relationship between the gradient and the depth is relatively stable, and the parameters of the depth function can be further optimized through calibration experiments.

[0049] It can be understood that the embodiments of the present invention do not strictly limit the filter size or operator parameters, and those skilled in the art can adjust the configuration according to the size of the organoid and the imaging conditions. For example, for organoids with a diameter less than 300 microns, the Sobel operator size can be appropriately reduced to improve the edge detection sensitivity. The obtained three-dimensional structure data of the organoid provides reliable support for drug sensitivity analysis, can accurately reflect the deformation characteristics, and meet the high-efficiency requirements of parallel processing of multiple samples.

[0050] S103. During the drug sensitivity experiment, aiming at the change characteristics of the volume and surface shape of the organoid, the optical flow method is used to track the position change of the key feature points on its surface and obtain the coordinates of the deformed feature points.

[0051] In the deformation tracking of the organoid in the drug sensitivity experiment, the dynamic changes of the surface feature points need to be concerned. For this purpose, first, the region growing method is used to extract the region with a curvature value greater than the preset threshold from the surface of the organoid as the sampling position of the feature points. The curvature value reflects the surface wrinkles and local morphological features. Usually, when the curvature exceeds 0.08, this region is regarded as a high-curvature region and is suitable as the tracking target. Taking an organoid with a diameter of 500 microns as an example, 200 to 300 feature points can be evenly selected, with a spacing of about 20 microns, to ensure coverage of the main morphological change regions. Subsequently, the Harris corner detector is used to screen the feature points, and the response threshold is set to 0.01 to ensure the stability and distinctiveness of the feature points, generating the initial surface feature point data.

[0052] S1031. In the feature point motion tracking stage, calculate the local gray gradient direction and amplitude according to the initial surface feature point data, and gradually reduce the resolution by performing four-layer pyramid decomposition on the collected image sequence to improve the calculation efficiency, where the base layer maintains the original resolution and the resolution is halved for each layer up; on this basis, apply the dense optical flow algorithm to calculate the motion vector field of the feature points, set the search window size to 16×16 pixels, which can cover the local deformation range of the organoid, and solve the displacement vector and velocity components by iteratively optimizing the optical flow constraint equation, generating the surface deformation feature data including the displacement field and velocity field, reflecting the change characteristics of the feature points in space and time.

[0053] For the surface deformation feature data, establish a Markov state transition matrix to describe the temporal dependence of the feature point positions. The state space of the matrix includes the position coordinates and velocity components, and the transition matrix is constructed by statistically analyzing the transition probabilities of the feature points between adjacent frames. To predict the morphological change trend of the organoid, use the Bayesian prediction algorithm to calculate the future position probability distribution in combination with historical data, set the prediction error threshold to 5 microns, and if the error exceeds this value, the prediction result is excluded to obtain reliable feature point position prediction data. This prediction method can capture the deformation trend in advance.

[0054] S1032. In the trajectory optimization stage, a particle filter is used to correct the predicted data of the feature point positions. The number of particles is set to 200, and each particle represents a possible state hypothesis. The state credibility is evaluated by calculating the particle weights, and the particles with weights lower than 0.01 are removed to reduce noise interference. Subsequently, the two-dimensional image coordinates are mapped to the three-dimensional space by using the three-dimensional orthogonal basis tensor projection. The projection matrix is generated from the camera calibration parameters, and the displacement of the surface grid nodes is calculated by triangulation to obtain the deformed spatial coordinate data of the feature points. The grid cells adopt a triangular structure to ensure the geometric consistency of the deformation calculation.

[0055] The coordinates of the deformed feature points are further optimized, and a regression neural network is used to refine the feature point trajectories. The network design includes 3 hidden layers with the number of neurons being 64, 32, and 16 respectively. The trajectory smoothness is optimized through multi-layer feature extraction and non-linear mapping. The deformation threshold is set to 30% of the original size. If the node displacement exceeds this value, it is regarded as abnormal and corrected by the regression model, which can effectively remove the abnormal displacements caused by light changes or culture medium disturbances and ensure the reliability of the coordinate data.

[0056] The application of the optical flow method in tracking benefits from its sensitivity to continuous motion. Dense optical flow assumes that the gray values are constant between adjacent frames, combines the smoothness constraint, and uses the variational method to solve the global optimal solution, which is more suitable for capturing the overall deformation of the organoid surface compared to sparse optical flow. Pyramid decomposition reduces the computational complexity through multi-scale analysis while retaining the detailed features. Experiments show that when the displacement of the feature points exceeds 50 microns, it usually corresponds to significant drug response deformations. At this time, the velocity field can further characterize the deformation rate and provide a dynamic basis for drug efficacy evaluation.

[0057] The advantage of combining the Markov matrix and Bayesian prediction lies in its ability to model time series data. The transition matrix captures the short-term dependencies of the feature points through statistical methods, while the Bayesian algorithm uses the prior distribution to predict the long-term trends. The two work together to improve the prediction accuracy. The introduction of the particle filter simulates multiple possible paths of the feature point trajectories through Monte Carlo sampling, and the weight update is based on the observation likelihood function, which can dynamically adapt to the non-linear characteristics of organoid deformation.

[0058] It can be understood that the embodiments of the present invention do not strictly limit the window size or the number of particles, and those skilled in the art can adjust the parameters according to the organoid size and experimental conditions. For example, for organoids with small deformations, the search window can be appropriately reduced to 12×12 pixels to improve the local accuracy. The finally obtained deformed feature point coordinate data completely records the morphological evolution process of the organoids under the action of drugs and provides high-precision input for three-dimensional structure analysis.

[0059] S104, matching the coordinates of the deformed feature points with the initial three-dimensional structure data of the organoid, and calculating the focal plane offset caused by the deformation of the organoid by least squares fitting.

[0060] First, the nearest neighbor matching algorithm is used to match the coordinates of the deformed feature points with the feature points in the initial three-dimensional structure data, and the point correspondence relationship is established by calculating the Euclidean distance between the two sets of feature points. Taking an organoid with a diameter of 400 microns as an example, the average distance between its surface feature points is about 20 microns. The rejection threshold is set to three times the average distance, that is, 60 microns. Matching points outside this range are considered abnormal and removed to generate reliable feature point matching data. Experimental verification shows that under normal culture conditions, about 90% of the feature point matching distance is less than 40 microns, which ensures the robustness of the matching results. Subsequently, the displacement of each feature point on the three-dimensional coordinate axis is calculated based on the matching data, and the weighted average method is used to determine the point deformation weight. The weight value reflects the significance of the local deformation, and the points with a weight greater than the preset threshold of 0.3 are selected as valid deformation data. Usually 60% to 70% of the points meet the conditions, mainly concentrated in the high curvature areas of the organoids. The deformation of these areas can better reflect the effect of the drug.

[0061] S1041. Based on the effective deformation data, a three-dimensional coordinate transformation matrix is ​​constructed to describe the spatial characteristics of the deformation. The matrix is ​​decomposed by the singular value decomposition method to extract the main transformation direction and amplitude, where the singular value represents the deformation contribution in each direction, and the first three largest singular values ​​usually account for more than 85% of the total change. The main transformation components are retained by setting a screening threshold, for example, components with singular values ​​less than 5% of the total are eliminated to obtain simplified spatial transformation feature data, thereby reducing noise interference and highlighting the dominant mode of deformation.

[0062] For the spatial transformation feature data, a five-layer deep neural network is designed to further explore the potential relationship between deformation and focal plane offset. The number of nodes in the network input layer is consistent with the feature dimension. For example, when a 150-dimensional feature vector is input, the number of nodes in the hidden layer decreases to 128, 64, 32, and 16 respectively, and the output layer generates the displacement components in the three directions of the focal plane. The network uses the ReLU activation function to enhance the nonlinear expression ability. The initial learning rate is set to 0.001 during training, and it decays to 0.8 times the original after every 1000 iterations to ensure convergence stability. After the focal plane displacement data is obtained through network prediction, principal component analysis is applied to extract the main offset direction, and the feature vector with a cumulative contribution rate of 90% is selected. Usually, the first 3 to 4 principal components can cover the main change trend and generate focal plane offset feature data. The key information is retained through dimensionality reduction while reducing the computational complexity.

[0063] A least - squares fitting equation system is established based on the focal - plane offset characteristic data, and the conjugate - gradient method is used to iteratively solve for the offset. The conjugate - gradient method gradually approaches the optimal solution by constructing an orthogonal direction sequence, and is particularly suitable for dealing with sparse matrix problems. The convergence threshold is set to 0.01 micrometers, and the iteration is terminated when the change in the solution vector for 5 consecutive iterations is less than this value. Experiments show that the focal - plane offset caused by deformation is most significant in the direction perpendicular to the culture plate, usually between 15 and 25 micrometers, while the offset parallel to the plane is smaller, about 5 to 10 micrometers, and the algorithm usually converges within 20 to 30 iterations. The accurate calculation of the offset can effectively correct the focal - plane position of the imaging system and avoid imaging blurring caused by the accumulation of deformation.

[0064] The weighted - average method considers the spatial distribution characteristics of the displacement vector in weight calculation. High - weight points usually correspond to significant change regions on the surface of the organoid under drug action, such as edge folds or volume expansion areas. This screening mechanism improves the representativeness of the data. Singular - value decomposition reveals the internal law of deformation through mathematical decomposition. Compared with directly calculating the displacement mean, it can better reflect the anisotropic characteristics of deformation. The introduction of a deep neural network further enhances the adaptability of the model to complex deformation patterns, captures the non - linear mapping relationship between displacement and focal - plane offset through multi - layer feature extraction, and ensures that the prediction results are close to the actual optical requirements.

[0065] It can be understood that the embodiments of the present invention do not make strict regulations on the number of network layers or the number of iterations. Those skilled in the art can adjust the parameters according to the complexity of organoid deformation. For example, for organoids with a larger deformation amplitude, the depth of the hidden layer can be increased to improve the fitting ability. The obtained focal - plane offset data provides a high - precision basis for subsequent imaging adjustment, effectively improving the imaging quality and drug - sensitivity analysis reliability in the scenario of parallel processing of multiple samples.

[0066] S105. According to the focal - plane offset, combined with the refractive index and absorption coefficient of the multi - well plate, use the ray - tracing algorithm to simulate the light - propagation path to correct the focal - plane position. At the same time, derive the phase - compensation matrix through wave - front aberration analysis to further optimize the focal - length adjustment accuracy of the imaging system.

[0067] Construct a three-dimensional grid structure of the multi-well plate based on the focal plane offset, and collect the refractive index and absorption coefficient data for each grid cell. Taking a 96-well plate as an example, the refractive index of the well positions usually ranges from 1.45 to 1.52, and the absorption coefficient fluctuates between 0.02 and 0.08. Calculate the distribution of material parameters between grid nodes by cubic spline interpolation method, and use 20 control points to generate a continuous optical property distribution. The interpolation accuracy is controlled within 0.1% by 5-fold cross-validation to ensure the spatial smoothness and authenticity of the parameters. For the obtained optical property distribution data of the multi-well plate, use the Monte Carlo ray tracing model to simulate the light propagation. Set 10,000 incident rays for each well position, and ensure the uniform light distribution by stratified sampling method. Backward tracing starts from the imaging plane, calculates the positions and incident angles of the intersection points of the rays and each interface of the multi-well plate, determines the refraction angle according to Snell's law, and records the total reflection event when the incident angle exceeds the critical angle of 41.8 degrees. Calculate the interface reflection loss in combination with the Franck-Condon principle, and the interfaces with reflectivity greater than 0.1 are marked as key interfaces, which are usually concentrated in the well wall area, and finally generate the light spatial propagation data and interface optical property data.

[0068] S1051. In the scattering and attenuation analysis, calculate the scattered light intensity distribution for the interface optical property data, construct a spatial model of scattering loss using a deep neural network. The network input includes position coordinates and local material parameters, and the output is the scattering coefficient. Optimize the model by dividing 80% of the training data and 20% of the validation data; the scattered intensity decays exponentially with depth, and the scattering increases significantly when the refractive index gradient exceeds 0.01 / μm. Based on this, establish the Lambert-Beer attenuation equation. The attenuation coefficient is about 0.05 / μm in the shallow layer of 0 to 200 μm and increases to 0.08 / μm in the deep layer of 200 to 500 μm. Calculate the light intensity compensation coefficient by iterative calculation and smooth it with a 50-μm Gaussian kernel to obtain the light intensity compensation data.

[0069] For the light intensity compensation data, construct an optimization equation set for the focal plane position and solve it using the conjugate gradient method. Set the initial step size to 1 μm and the convergence threshold to 0.1 μm. Convergence occurs after 15 to 20 iterations. Eliminate abnormal solutions by standard deviation analysis, and control the focal plane position error within 2 μm. The advantage of Monte Carlo ray tracing is that it can accurately reflect the influence of the material inhomogeneity of the multi-well plate on light propagation by simulating a large number of random samples of light paths, and is more suitable for complex optical environments compared with deterministic methods. Experiments show that the corrected focal plane position effectively reduces the imaging distortion caused by material distribution differences, and the focal plane drift is reduced by about 70%, significantly improving the stability of long-term observation.

[0070] S1052. In wavefront aberration processing, a double-beam interferometer is used to measure the refractive index distribution of the porous plate grid cells with a wavelength of 633 nm. The optical path change is calculated through the phase difference of the interference fringes. The 128-point fast Fourier transform is applied to extract the frequency components. The cut-off frequency is set to 1 / 4 of the sampling frequency to filter out high-frequency noise, and a refractive index distribution map is generated. The wavefront aberration is decomposed using third-order Zernike polynomials, which include 9 basis functions to characterize the aberration types. The coefficients are calculated by the least squares method, and noise terms are removed with a screening threshold of 0.05. The coefficient of the defocus term is usually between 0.2 and 0.4, and the astigmatism term is between 0.1 and 0.2, obtaining the wavefront aberration feature vector.

[0071] Based on the wavefront aberration feature vector, a three-layer neural network is constructed to optimize the scattering parameters. The 9 nodes in the input layer correspond to the aberration coefficients, the hidden layer has 16 nodes, and the 4 nodes in the output layer generate parameters such as scattering angles and intensities. The hyperbolic tangent activation function is used. The training uses the stochastic gradient descent method with a batch size of 32, and the initial learning rate is 0.01, which decays to 0.8 times every 1000 iterations. Combining with the path random distribution model, the zero-order Bessel function is used to describe the radial distribution of the scattering intensity, and the fourth-order high-order term is calculated through the recurrence formula. In the small-angle scattering region, the zero-order term dominates, and the attenuation rate is about 2 times that of the cosine function. The scattering distribution is convolved with the absorption coefficient to generate phase compensation data. The singular value decomposition is used to extract the phase compensation matrix, and the first 5 singular value components are retained, with a contribution rate of over 95%. Based on this matrix, the focal plane offset correction parameters are updated and optimized by the conjugate gradient method, with a gradient threshold of 0.01 and a maximum of 50 iterations, and the final error is controlled within 2 microns.

[0072] Wavefront aberration analysis provides a theoretical basis for phase compensation by quantifying the perturbation of the optical path by the refractive index distribution. The decomposition of Zernike polynomials reveals the main aberration types caused by the inhomogeneity of the porous plate material, and the neural network enhances the prediction accuracy of the scattering parameters through non-linear mapping. Experimental verification shows that the wavefront aberration amount is reduced by 75% after phase compensation, significantly improving the real-time performance of focal plane tracking and the imaging clarity.

[0073] It can be understood that the embodiments of the present invention do not impose fixed restrictions on the number of light rays or the grid size. Those skilled in the art can adjust the parameters according to the specifications of the porous plate. For example, for a 384-well plate with a smaller pore diameter, the light ray density can be appropriately increased to improve the simulation accuracy. The focal plane position correction data provides reliable support for imaging optimization, ensuring the accurate capture of the morphological characteristics of the organoids in the drug sensitivity experiment.

[0074] S106. Based on the corrected focal plane position, a fast focusing algorithm is used to adjust the focal length of the optical imaging system to re-acquire the organoid images, and high-precision focusing is achieved through beam propagation characteristic analysis and multi-level optimization, ensuring that the imaging quality meets the requirements of the drug sensitivity experiment.

[0075] Construct the Gaussian beam propagation equation according to the corrected focal plane position to describe the spatial distribution characteristics of the beam, and generate the focal length adjustment curve by calculating the mapping relationship between the beam cross-sectional radius and the optical path length. The beam cross-sectional radius reaches the minimum value at the focus, usually about 20 microns, and shows a hyperbolic variation law along the propagation direction. To optimize the curve accuracy, the gradient descent method with an adaptive learning rate is used for iteration. The initial learning rate is set to 0.01. If the gradient change in 5 consecutive iterations is less than 0.001, it is reduced to 0.8 times to obtain the initial focal length adjustment parameter. Then, use the Michelson interferometer to measure the optical path difference. The optical path difference between the reference arm and the measurement arm is controlled within 200 microns. Extract the numerical values through the phase distribution of the interference fringes, and then perform multi-scale decomposition using the discrete wavelet transform. Select the db4 wavelet basis function and decompose it to 4 layers. Select the reconstruction coefficients of the second and third layers because their energy ratio exceeds 85% to generate the optical path compensation data.

[0076] S1061. In the focal parameter optimization stage, input the optical path compensation data into a five-layer convolutional neural network. The network structure is designed with 32 nodes in the input layer, 64, 128, and 64 nodes in the three hidden layers respectively, and 32 nodes in the output layer. Residual connections are used to span adjacent layers to alleviate the problem of gradient disappearance; the network outputs the focal parameter curve, which is numerically normalized to the interval from -1 to 1 after batch normalization, facilitating the adjustment of the actuator, and at the same time enhancing the adaptability of the algorithm to different optical path changes.

[0077] For the standardized focal curve, divide the focal length adjustment range into 16 sub-intervals by the uniform division method, construct a state transition matrix to describe the adjustment cost between adjacent intervals, and comprehensively consider the power consumption and adjustment time of the piezoelectric ceramic actuator. Use the forward dynamic programming algorithm to calculate the minimum cost path, recursively deduce from the initial state to the final state layer by layer, and generate the focusing adjustment sequence. Based on this sequence, control the piezoelectric ceramic actuator to adjust the lens group spacing. The displacement resolution of the actuator reaches 0.1 micron, the maximum stroke is 500 microns, and the response time is less than 1 millisecond, ensuring the real-time and accuracy of the adjustment.

[0078] Use the Gaussian kernel support vector machine to evaluate the clarity of the acquired images. The kernel parameter σ is set to 2.0, and the penalty factor C is optimized to 10 through cross-validation. Calculate the variance and gradient features of the images, where the variance represents the contrast and the gradient reflects the edge sharpness. Construct a focal length fine-tuning function based on the clarity evaluation data, calculate the focal length deviation by the weighted least squares method, and perform closed-loop adjustment in combination with a proportional-integral controller. Set the proportional coefficient Kp to 0.8, the integral coefficient Ki to 0.2, and the dead zone threshold to 0.5 micron to avoid system oscillation and improve the steady-state accuracy. Experimental verification shows that when the focal length deviation is less than the dead zone threshold for 5 consecutive times, the system reaches stability, and the entire focusing process is usually completed within 200 milliseconds, meeting the high-efficiency requirements of real-time imaging.

[0079] The Gaussian beam propagation equation provides a theoretical basis for focal length adjustment by analyzing the relationship between the beam waist position and the propagation distance, and its accuracy directly affects the subsequent imaging effect. The multi-scale decomposition of the discrete wavelet transform effectively separates the low-frequency trend and high-frequency noise in the optical path difference, and can retain key features better than traditional filtering. The residual structure of the convolutional neural network enhances the deep feature extraction ability, enabling the focusing parameters to dynamically adapt to the complex changes in the deformation of the organoid and the optical properties of the multi-well plate. The dynamic programming algorithm reduces the energy consumption and time cost during the adjustment process through global optimization, ensuring the focusing efficiency. The clarity evaluation model of the support vector machine uses non-linear mapping to convert the image features into quantifiable clarity indicators, providing data-driven support for fine-tuning the focal length.

[0080] It can be understood that the embodiments of the present invention do not strictly limit the number of network layers or the number of interval partitions, and those skilled in the art can adjust the parameters according to the resolution of the imaging system and the size of the organoid. For example, for a high-resolution system, the decomposition level can be increased to improve the optical path compensation accuracy. The adjusted focal length significantly improves the clarity of the re-acquired organoid images, providing high-quality image data support for three-dimensional reconstruction and drug sensitivity analysis.

[0081] S107. Use a three-dimensional reconstruction algorithm stacked layer by layer to register and splice the re-acquired image sequence to generate the complete three-dimensional structure data after the deformation of the organoid.

[0082] Construct a corner response function based on the re-acquired image sequence to extract feature points, calculate the corner response value using the Harris operator, and set the threshold to 0.01 to filter out stable corners. Taking an organoid with a diameter of 500 microns as an example, 200 to 300 feature points can be extracted, and redundant points are removed by non-maximum suppression with a radius of 5 pixels to ensure uniform distribution. Then, calculate the local binary pattern features using a circular sampling window, set the window radius to 8 pixels, sample 16 points to generate a gray-scale binary code, and form the image feature point data. This method enhances the distinguishability of feature points through local texture information, providing a reliable basis for subsequent registration. For the feature point data, construct a cross-correlation response function, calculate the translation amount and rotation angle between adjacent images by the phase correlation method, combine the fast Fourier transform to convert the spatial domain data into the frequency domain, set the transformation resolution to 256×256 points, the peak position reflects the translation amount, the peak amplitude evaluates the registration reliability, and when the normalized cross-correlation coefficient is greater than 0.8, the registration is considered effective, obtaining the image registration data. Experiments show that the translation amount between adjacent layers is usually less than 20 microns, and the rotation angle does not exceed 5 degrees, fully reflecting the continuity of the organoid deformation.

[0083] S1071. In the interlayer registration stage, calculate the correspondence of feature points based on the image registration data, use the nearest neighbor distance ratio method to screen the matching pairs, set the distance threshold to 0.7 to eliminate the wrong matches, optimize the interlayer transformation matrix by the least squares method, decompose the matrix into translation, rotation, and scaling components, use singular value decomposition to extract the main transformation modes, retain the first three singular value components, whose contribution rate usually exceeds 95%, and generate accurate spatial transformation parameters.

[0084] For the spatial transformation parameters, construct a five-layer deep neural network to extract image features. The convolution kernel size is uniformly 3×3, and the number of channels is 32, 64, 128, 64, and 32 in sequence. Through a four-layer spatial pyramid structure, increase the pooling kernel size from 2×2 to 16×16 to capture multi-scale spatial information. Adopt skip connections to fuse shallow details and deep semantic features to generate feature fusion data and enhance the model's perception ability of deformation details. Based on the feature fusion data, perform deconvolution reconstruction, use bilinear interpolation to restore the image resolution, avoid the checkerboard effect, and then apply the Canny edge detection operator to extract the structural contour. Set the high and low threshold ratio to 2:1 to ensure the continuity of the contour. The registered contour points form three-dimensional contour data.

[0085] Based on the three-dimensional contour data, construct hexahedral mesh elements with an element size of about 10 microns. Optimize the node distribution through the Markov energy function. The energy function combines the data term and the smoothing term. The data term measures the fitting degree of the nodes to the observed data, and the smoothing term constrains the continuity and uniformity of the mesh form. Use the Delaunay criterion for triangulation to ensure that the minimum angle of the triangle is greater than 30 degrees to avoid the influence of long and narrow triangles on the mesh quality. Combine the Laplace operator for surface smoothing. Adjust the node positions towards the local average distribution in each iteration, and control the number of iterations within 10 times. Finally, generate the complete three-dimensional structure data of the organoid with a spatial resolution reaching the sub-micron level.

[0086] The application of the phase correlation method in registration benefits from its robustness to brightness changes and noise. It achieves sub-pixel accuracy through frequency domain peak localization, significantly improving the interlayer alignment efficiency. The multi-scale feature fusion of the deep neural network retains information at different levels through skip connections and is more adaptable to the irregularity of organoid deformation compared to single-scale analysis. The Markov energy optimization balances the data fidelity and the structural smoothness through global energy minimization, ensuring that the reconstructed model is both accurate and natural. Experimental verification shows that this method can effectively capture the spatial details of the organoid when dealing with complex deformations, such as surface wrinkles or volume changes.

[0087] It is understandable that the embodiments of the present invention do not impose strict regulations on the number of feature points or the number of network layers, and those skilled in the art can adjust the parameters according to the image resolution and the complexity of the organoid. For example, for a high-density image sequence, the number of pyramid layers can be increased to enhance the feature extraction ability. The obtained three-dimensional structure data completely reflects the morphological evolution of the organoid in the drug sensitivity experiment.

[0088] S108. Extract the volume, surface area, and cell distribution density features of the organoid based on the complete three-dimensional structure data, and generate a morphological feature report for the drug sensitivity experiment through multi-level calculation and analysis.

[0089] According to the vertex coordinates and patch information of the three-dimensional structure data, a tetrahedral meshing function is constructed using the least squares criterion to calculate the volume of the organoid. The meshing size threshold is set to 20 microns to balance accuracy and computational efficiency. The voxel volume of each tetrahedral element is calculated through the three-point Gaussian quadrature formula, and the integration point weights are 0.25, 0.25, and 0.5 respectively, with the error controlled within 0.1%. Taking an organoid with a diameter of 500 microns as an example, the grid is usually divided into 8000 to 10000 tetrahedral elements, and the initial volume data is finally obtained. For the grid surface, the barycentric coordinate mapping method is used to calculate the area of the triangular patches. The geometric properties of each patch are determined by the vertex coordinates, and then a surface weighting function is constructed by combining the mean curvature and the Gaussian curvature. The area calculation deviation is adjusted by multiplying the area of the region where the mean curvature is greater than 0.1 or the Gaussian curvature is less than -0.01 by a correction factor of 1.1 to 1.3 to generate the initial surface area data. Then, a five-layer convolutional neural network is constructed to process the surface area data. The network structure includes an input layer with 32 channels, three hidden layers with 64, 128, and 64 channels respectively, and an output layer with 32 channels. A batch normalization layer and a ReLU activation function are connected after each layer. The feature dimension is compressed through a pooling layer and multi-scale features are fused using a fully connected layer to obtain the surface area feature data. This neural network design can effectively extract the deep features of the surface morphology and improve the robustness of feature expression.

[0090] For the three-dimensional structural grayscale atlas, the Otsu method is used to calculate the adaptive segmentation threshold to separate the nucleus region. This method automatically determines the optimal threshold by analyzing the between-class difference of the grayscale histogram, adapting to the image characteristics under different lighting conditions. The twenty-six neighborhood criterion is used to label the three-dimensional connected regions, and the complete spatial structure of the nucleus is identified through a recursive scanning algorithm. The centroid position of each connected region is calculated and the quantity is counted to generate the nucleus position data. Taking an organoid with a volume of 0.1 cubic millimeters as an example, it usually contains 1,000 to 2,000 nuclei. Based on the nucleus position data, a spatial density function is constructed, and the kernel density estimation method is used to calculate the local density distribution. The Gaussian kernel function is selected, and the bandwidth is set to 30 microns, which not only ensures a smooth distribution but also retains local details. Subsequently, the density gradient field is calculated through radial basis support vector regression, the kernel parameter σ is set to 50 microns, the penalty factor C is 10, and the cross-validation error is less than 5%, obtaining the cell density characteristic data. The local maximum value in the gradient field reflects the spatial characteristics of the cell aggregation region.

[0091] S1081. In the feature integration stage, the initial volume data, surface area feature data, and cell density feature data are standardized. The top five morphological parameters with the highest feature contribution rates are extracted through principal component analysis, usually including the volume change rate, surface fold degree, cell density uniformity, morphological complexity, and growth activity, and the cumulative contribution rate reaches more than 90%. The cubic polynomial fitting method is used to construct the feature description function, and the least squares method is used to optimize the fitting parameters. The root mean square error of the residuals is controlled within 3% to generate the morphological feature report data, including indicators such as the volume size and change trend, the surface area to volume ratio, the local curvature distribution, the spatial characteristics and aggregation degree of cell density.

[0092] The tetrahedral meshing optimizes the mesh division through the least squares criterion to ensure the geometric consistency of the voxel units. The three-point Gaussian integration accurately calculates the volume through numerical integration, and is more adaptable to the irregular shape of the organoid compared to simple voxel summation. The curvature weighting function in the surface area calculation corrects the area deviation under the flat assumption by identifying the folded regions, making the result closer to the true surface characteristics. The multi-scale feature fusion of the neural network enhances the perception ability of surface details through pooling and fully connected operations, especially when significant folds occur on the surface of the organoid under the action of drugs. The combination of kernel density estimation and support vector regression in cell density analysis reveals the dynamic changes in cell distribution through non-parametric estimation and regression modeling, providing a spatial basis for evaluating the impact of drugs on cell proliferation or apoptosis.

[0093] It can be understood that the embodiments of the present invention do not impose fixed restrictions on the dissection size or the number of network channels. Those skilled in the art can adjust the parameters according to the size of the organoids and the experimental requirements. For example, for small organoids, the dissection threshold can be reduced to improve the accuracy. The generated morphological feature report quantifies the deformation characteristics of the organoids through multi-dimensional indicators, providing a comprehensive and reliable analysis basis for the drug sensitivity experiment, and significantly improving the evaluation efficiency and scientificity of the drug response.

[0094] The present invention provides an automatic organoid drug sensitivity detection system for parallel processing of multiple samples, mainly including:

[0095] An optical path length calculation module, configured to obtain the optical path length data of each well position of the cell culture multi-well plate, and calculate the initial focal plane position of each well position through a preset optical model;

[0096] An initial three-dimensional structure reconstruction module, configured to collect the organoid images of each well position according to the initial focal plane position, extract the organoid contour by using an edge detection algorithm, and map the two-dimensional contour to a three-dimensional space model by using a shape matching method to obtain the initial three-dimensional structure data of the organoid;

[0097] A feature point tracking module, configured to, during the drug sensitivity experiment, track the position changes of the key feature points on the surface of the organoid by using an optical flow method according to the volume and surface shape change characteristics of the organoid, and obtain the coordinates of the deformed feature points;

[0098] A focal plane offset calculation module, configured to match the coordinates of the deformed feature points with the initial three-dimensional structure data of the organoid, and calculate the focal plane offset caused by the deformation of the organoid by using the least squares fitting method;

[0099] A ray tracing simulation module, configured to, according to the focal plane offset, combine the refractive index and absorption coefficient of the multi-well plate, and simulate the refraction and scattering paths of light in the multi-well plate by using a ray tracing algorithm to obtain the corrected focal plane position;

[0100] A focal length adjustment module, configured to adjust the focal length of the optical imaging system according to the corrected focal plane position by using a fast focusing algorithm, and re-collect the organoid images;

[0101] A three-dimensional reconstruction module, configured to register and splice the re-collected image sequences by using a three-dimensional reconstruction algorithm of layer-by-layer stacking to obtain the complete three-dimensional structure data of the deformed organoid;

[0102] A morphological feature generation module, configured to extract the volume, surface area, and cell distribution density features of the organoid according to the complete three-dimensional structure data, and generate a morphological feature report for the drug sensitivity experiment.

[0103] It should be further noted that, for the various specific technical features described in the above specific embodiments, they can be combined in any appropriate manner without contradiction. To avoid unnecessary repetition, the present invention will not separately describe various possible combination methods. In addition, any combination can be made between various different embodiments of the present invention, as long as it does not violate the idea of the present invention, it should also be regarded as the content disclosed by the present invention.

Claims

1. An automatic detection method for organoid drug sensitivity with parallel processing of multiple samples, characterized in that, The method includes: Obtaining the optical path length data of each well position of a cell culture multi-well plate and determining the initial focal plane position of each well position; According to the initial focal plane position, acquiring the organoid images of each well position, extracting the organoid contour using an edge detection algorithm, and mapping the two-dimensional contour to a three-dimensional space model using a shape matching method to obtain the initial three-dimensional structure data of the organoid; During the drug sensitivity experiment, aiming at the characteristics of the volume and surface shape change of the organoid, using the optical flow method to track the position change of the key feature points on the organoid surface to obtain the coordinates of the deformed feature points; Matching the coordinates of the deformed feature points with the initial three-dimensional structure data of the organoid, and calculating the focal plane offset caused by the organoid deformation through least squares fitting; According to the focal plane offset, combining the refractive index and absorption coefficient of the multi-well plate, using the ray tracing algorithm to simulate the refraction and scattering paths of light in the multi-well plate to obtain the corrected focal plane position; Adopting a fast focusing algorithm, adjusting the focal length of the optical imaging system according to the corrected focal plane position, and re-acquiring the organoid images; Using a three-dimensional reconstruction algorithm of layer-by-layer stacking to register and splice the re-acquired image sequence to obtain the complete three-dimensional structure data of the deformed organoid; According to the complete three-dimensional structure data, extracting the volume, surface area, and cell distribution density characteristics of the organoid, and generating a morphological feature report of the drug sensitivity experiment.

2. The method according to claim 1, wherein The obtaining the optical path length data of each well position of a cell culture multi-well plate and determining the initial focal plane position of each well position includes: Collecting the optical path length data of the well positions of the multi-well plate through an optical sensor, and screening the optical path length data according to a preset upper limit value and lower limit value of light intensity to obtain an initial well position optical data set; Using an optical path sensor to scan the multi-well plate layer by layer, calculating the variance of the gray value of each layer according to the scanned data to obtain a clarity index, and selecting the layer corresponding to the maximum clarity index as the focal plane candidate layer as the initial focal plane position of each well position.

3. The method according to claim 1, wherein The collecting the organoid images of each well position according to the initial focal plane position, extracting the organoid contour using an edge detection algorithm; and mapping the two-dimensional contour to a three-dimensional space model using a shape matching method to obtain the initial three-dimensional structure data of the organoid includes: Adopting a multi-well plate scanning method for the organoid focal plane position data to obtain an image sequence, and processing the image sequence through Gaussian filtering and Laplace operator to obtain enhanced organoid image data; Extracting a contour point set using a multi-scale sobel operator according to the enhanced organoid image data, and performing double-threshold connection processing on the contour point set to obtain edge contour feature data; Performing polar coordinate transformation and depth estimation on the edge contour feature data, and using the Delaunay triangulation algorithm to perform grid processing on the contour feature points to obtain spatial grid structure data; Applying a grid smoothing algorithm based on Laplace coordinates and an iterative grid deformation algorithm to the spatial grid structure data, and calculating the vertex displacement of the grid according to the least squares method to obtain the three-dimensional structure data of the organoid.

4. The method according to claim 1, wherein During the drug sensitivity experiment, according to the volume and surface shape change characteristics of the organoids, the optical flow method is used to track the position changes of the key feature points on the surface of the organoids, and the coordinates of the deformed feature points are obtained, including: The region where the surface curvature value of the organoids is greater than the preset curvature threshold is obtained by the region growing method, and the Harris corner detector with a response value greater than the preset response threshold is used to screen the feature points in the region to obtain the initial surface feature point data; The local gray gradient direction and gradient amplitude are calculated according to the initial surface feature point data, and the dense optical flow algorithm with a window size of a fixed pixel value is used to calculate the feature point motion vector field to obtain the surface deformation feature data; A Markov state transition matrix is established for the surface deformation feature data, and the Bayesian prediction algorithm is used to predict the morphological change trend of the organoids. If the prediction error is less than the preset error threshold, the feature point position prediction data is obtained; A fixed number of particle filters are used to correct the trajectory of the feature point position prediction data, and the spatial coordinates of the feature points are calculated by three-dimensional orthogonal basis tensor projection to obtain the coordinates of the deformed feature points.

5. The method according to claim 1, wherein The coordinates of the deformed feature points are matched with the initial three-dimensional structure data of the organoids, and the focal plane offset caused by the deformation of the organoids is calculated by the least squares fitting, including: The nearest neighbor matching algorithm is used to pair the coordinates of the deformed feature points with the initial three-dimensional structure feature points, and the point position correspondence relationship is established according to the Euclidean distance between the feature points. The points with a matching distance greater than the average distance threshold between points are removed to obtain the feature point matching data; The displacement of the point position in the three coordinate axes is calculated according to the feature point matching data, and the weighted average method is used to calculate the point position deformation weight. The points with a weight less than the preset weight value are screened to obtain the effective deformation data; A three-dimensional coordinate transformation matrix is constructed for the effective deformation data, and the main transformation direction and transformation amplitude are extracted by singular value decomposition. The main transformation components are retained by setting the singular value screening threshold to obtain the spatial transformation feature data; A five-layer deep neural network is constructed according to the spatial transformation feature data, and the principal component analysis is used to select the eigenvectors with a cumulative contribution rate reaching the preset threshold as the main offset direction to obtain the focal plane offset.

6. The method according to claim 1, wherein According to the focal plane offset, combined with the refractive index and absorption coefficient of the multi-well plate, the ray tracing algorithm is used to simulate the refraction and scattering paths of light in the multi-well plate to obtain the corrected focal plane position, including: A three-dimensional grid structure of the multi-well plate is constructed according to the focal plane offset, and the cubic spline interpolation method is used to calculate the material parameter distribution between the grid nodes to obtain the optical property distribution data of the multi-well plate; The stratified sampling method is used to generate the spatial distribution of incident light for the optical property distribution data of the multi-well plate, and the propagation direction of light in the multi-well plate is calculated by backtracking to obtain the light spatial propagation data; The Snell's law is used to calculate the interface refraction angle according to the light spatial propagation data, and the Frank-Condon principle is used to calculate the interface reflection loss to obtain the interface optical property data; Establish a Lambert-Beer attenuation equation for the interface optical characteristic data, and solve the optimization equation set of the focal plane position by the conjugate gradient method to obtain the corrected focal plane position; it also includes: obtaining the refractive index distribution data of the porous plate, calculating the wavefront distortion amount of the transmitted light, adjusting the scattering parameters of the preset path random distribution model, fusing the path random distribution model with the absorption coefficient, deriving the phase compensation data, and updating the correction parameters of the focal plane offset, specifically including: obtaining the refractive index data of the porous plate structure, extracting the frequency components of the refractive index data by fast Fourier transform, and performing high-frequency filtering according to the preset cut-off frequency to obtain the refractive index distribution map.

7. The method according to claim 1, characterized in that The fast focusing algorithm is adopted to adjust the focal length of the optical imaging system according to the corrected focal plane position, and re-collect the organoid images, including: Calculate the mapping relationship between the beam cross-sectional radius and the optical path length according to the Gaussian beam propagation equation to obtain the focal length adjustment curve; Use a Michelson interferometer to measure the optical path difference value corresponding to the focal length adjustment curve, and perform multi-scale decomposition on the optical path difference value by discrete wavelet transform to obtain the optical path compensation data; Input the optical path compensation data into a convolutional neural network, and the output layer of the convolutional neural network generates a focusing parameter curve, and obtain a normalized focusing curve through batch normalization processing; After controlling the piezoelectric ceramic actuator to adjust the lens group spacing according to the normalized focusing curve, re-collect the organoid images.

8. The method according to claim 1, wherein The three-dimensional reconstruction algorithm using layer-by-layer stacking is used to register and splice the re-collected image sequence to obtain the complete three-dimensional structure data of the deformed organoid, including: Construct a corner response function according to the image sequence, and extract feature point coordinate data by the corner response function for a preset corner response threshold; Based on the feature point coordinate data, calculate the translation amount and rotation angle, and perform spectral transformation on the translation amount and rotation angle by fast Fourier transform to obtain the image registration data; Calculate the corresponding relationship of feature points according to the image registration data, establish a matching distance threshold through the corresponding relationship of feature points, and optimize the matching distance threshold by the least square method to obtain the interlayer transformation matrix; Extract image features for the interlayer transformation matrix, and perform multi-scale fusion on the image features by skip connection to obtain feature fusion data; Extract the structure contour according to the feature fusion data, and perform spatial registration on the contour points to obtain the complete three-dimensional structure data of the deformed organoid.

9. The method according to claim 1, wherein According to the complete three-dimensional structure data, extract the volume, surface area and cell distribution density characteristics of the organoid, and generate a morphological feature report for the drug sensitivity experiment, including: According to the three-dimensional structure vertex coordinates and patch information, construct a tetrahedral meshing function using the least square criterion, and calculate the voxel unit volume by the three-point Gaussian integration formula to obtain the initial volume data; For the initial volume data, calculate the triangular patch area by barycentric coordinate mapping to obtain the surface area data; According to the surface area data, process the feature data through a batch normalization layer, compress the feature dimension by a pooling layer and fuse multi-scale features by a fully connected layer to obtain the surface area feature data; For the surface area feature data and the initial volume data, a feature description function is constructed by cubic polynomial fitting to obtain morphological feature report data.

10. An organoid drug sensitivity automatic detection system for parallel processing of multiple samples, characterized in that, The system includes: An optical path length calculation module, configured to obtain optical path length data of each well position of a cell culture multi-well plate, and calculate the initial focal plane position of each well position through a preset optical model; An initial three-dimensional structure reconstruction module, configured to collect organoid images of each well position according to the initial focal plane position, extract the organoid contour using an edge detection algorithm, and map the two-dimensional contour to a three-dimensional space model using a shape matching method to obtain the initial three-dimensional structure data of the organoid; A feature point tracking module, configured to, during the drug sensitivity experiment, track the position change of key feature points on the surface of the organoid using an optical flow method according to the volume and surface shape change characteristics of the organoid to obtain the coordinates of the deformed feature points; A focal plane offset calculation module, configured to match the coordinates of the deformed feature points with the initial three-dimensional structure data of the organoid, and calculate the focal plane offset caused by the deformation of the organoid through least squares fitting; A ray tracing simulation module, configured to, according to the focal plane offset, combine the refractive index and absorption coefficient of the multi-well plate, and simulate the refraction and scattering paths of light in the multi-well plate using a ray tracing algorithm to obtain the corrected focal plane position; A focal length adjustment module, configured to use a fast focusing algorithm to adjust the focal length of the optical imaging system according to the corrected focal plane position and re-collect organoid images; A three-dimensional reconstruction module, configured to use a three-dimensional reconstruction algorithm of layer-by-layer stacking to register and splice the re-collected image sequence to obtain the complete three-dimensional structure data of the deformed organoid; A morphological feature generation module, configured to extract the volume, surface area, and cell distribution density features of the organoid according to the complete three-dimensional structure data and generate a morphological feature report of the drug sensitivity experiment.

Images