A method for measuring the direction of movement of flagellate algae under flow conditions

Abstract

The application discloses a kind of methods for measuring the movement direction of flagellum algae under flow conditions, the method comprises the following steps: step 1, image sequence acquisition and processing;Step 2, cell image cropping alignment;Step 3, cell image pre-processing;Step 4, create fan ring detection area;Step 5, construct timing frame;Step 6, time sequence difference feature calculation;Step 7, confidence calculation and movement direction determination.The method for measuring the movement direction of flagellum algae under flow conditions disclosed in the application overcomes the technical problems that cannot be identified due to poor imaging quality of flagellum and loss of instantaneous information in moving water environment, and realizes high-precision, automated and continuous measurement of the movement direction of flagellum algae through time sequence dynamic feature analysis and biological motion characteristic determination.

Description

Technical Field

[0001] This invention belongs to the field of environmental hydrodynamics technology, and in particular relates to a method for measuring the direction of movement of flagellates under flowing conditions. Background Technology

[0002] Dyazoa blooms are a typical algal bloom phenomenon widely occurring in riverine reservoirs, estuaries, bays, and other water bodies. Dyazoa can actively swim by the beating of their flagella. This autonomous movement reflects the cell's response strategy to external stimuli (such as light, chemical gradients, gravity, and flow fields) and its internal motion regulation mechanism, directly determining its migration, aggregation characteristics, and spatiotemporal distribution patterns in water bodies. Therefore, accurately measuring the movement direction of dyazoa is crucial for revealing the dynamic mechanisms of bloom formation and dissipation and for achieving effective prediction and forecasting.

[0003] Shear flow is prevalent in both natural and artificial aquatic environments. In this context, algal cell movement is a complex result of the combined effects of their own dynamic drive, flow-driven motion, and the torque exerted by the environmental flow field. The shear flow field significantly influences cell orientation and rotation, causing their movement to exhibit dynamic and complex kinetic characteristics. Therefore, developing a technology capable of automatically and precisely analyzing the movement direction of flagellates under flowing conditions is crucial for revealing the mechanisms of water-algae interactions.

[0004] The movement of flagellates primarily relies on the propulsion of their cell flagella. Their cells are typically ellipsoidal, and the direction of movement is parallel to their long axis. However, this long axis has two opposing directions: one from the centroid to the location of the flagella, and the other from the centroid to the opposite direction. Therefore, the actual direction of cell movement must be determined based on the spatial position of the flagella relative to the cell's long axis, specifically depending on the flagellar propulsion pattern: for algal cells that rely on the pulling force generated by the flagella, the direction of movement is the same as the direction of the pulling force, i.e., the cell moves along the centroid towards the flagella; for algal cells that rely on the pushing force generated by the flagella, the direction of movement is the same as the direction of the pushing force, i.e., the cell moves along the centroid towards the opposite direction. Thus, accurately determining the spatial position of the flagella relative to the cell body is fundamental to determining its true direction of movement.

[0005] However, achieving automatic identification of flagella position faces significant technical challenges. Algal cell size is typically on the micrometer scale, while flagella diameter is only on the nanometer scale, requiring high-magnification imaging to observe both the cell body and flagella. Although current image recognition technologies can identify the cell body and its long axis direction relatively accurately, direct and automated morphological identification remains extremely challenging due to the tiny size of flagella, low contrast with the background, and susceptibility to motion blur during flow. This constitutes a fundamental technical bottleneck in motion direction detection. Furthermore, under the influence of shear flow fields, flagellated algal cells undergo rapid rotation. This rotation leads to two key problems: first, the constantly changing cell posture makes it easier for the flagella to detach from the microscope's focal plane during high-speed movement, resulting in blurred, low-contrast, or even completely invisible flagella morphology in single or multiple frames; second, the angular velocity of cell rotation is significantly accelerated in the flow field, further exacerbating the difficulty and instability of flagella imaging.

[0006] Traditional morphological recognition methods based on single-frame images rely entirely on the clear visibility of flagella in each frame. If flagella lose information due to detachment from the focal plane in one or more frames, the method cannot determine the cell's direction of motion at that moment, resulting in data interruption and loss of valuable information. This inability to make judgments due to momentary information loss is particularly pronounced in dynamic, shear-effect-laden flow environments, representing another limitation that current technologies struggle to overcome.

[0007] Currently, the detection of the movement direction of dinoflagellates mainly relies on two types of traditional methods:

[0008] Trajectory-velocity direction method: This method approximates the direction of motion by tracking the cell's trajectory and approximating the direction of velocity of the centroid. However, this method is only applicable to still or weakly flowing environments. In the presence of a significant background flow field (especially shear flow), cell motion is a superposition of its own swimming and environmental transport, and the velocity direction can no longer represent its own driving direction, causing the method to fail.

[0009] Manual observation and judgment method: This method involves researchers manually observing the long axis direction of cells and the position of flagella using microscopic imaging. This method is highly subjective, inefficient, and cannot perform large-scale statistical analysis, making it unsuitable for precise and reproducible studies on large numbers of dynamically changing samples.

[0010] Therefore, under flow conditions, especially in flow fields with significant shear, how to achieve automatic and accurate measurement of the direction of movement of dinoflagellates has become a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0011] The purpose of this invention is to provide a method for measuring the direction of movement of dinoflagellates under flowing conditions, so as to solve the above-mentioned technical problems.

[0012] To achieve the above objectives, the present invention provides the following technical solution:

[0013] This invention discloses a method for measuring the direction of movement of dinoflagellates under flowing conditions, the method comprising the following steps:

[0014] Step 1: Image sequence acquisition and processing: Image sequences recording the effective movement trajectories of flagellate cells are acquired using the experimental setup. The image sequences are then processed to identify the cell bodies of the algae, fit them into ellipses, and save the cell ellipse parameters, including the centroid coordinates, the length of the major and minor axes, and the direction of the major axis.

[0015] Step 2, Cell Image Cropping and Alignment: Read the cell ellipse parameters, set a standard image size that can encompass the entire cell, and translate the centroid of each cell to the center of the image to achieve spatial alignment of the cell image;

[0016] Step 3, Cell Image Preprocessing: The cropped and aligned cell image is preprocessed, including Gaussian blur noise reduction, background estimation and removal, and contrast enhancement.

[0017] Step 4: Create fan-ring detection regions: For the preprocessed cell image, fan-ring detection regions centered on the cell centroid are created at both ends of the cell's long axis. Specifically, the fan-ring detection regions are: the inner radius avoids the cell body, the outer radius covers the flagellar activity area, and symmetrical regions extend to both sides at a certain angle with the cell's long axis as the center; thus, a front-end fan-ring detection region and a rear-end fan-ring detection region are formed at both ends of the cell's long axis.

[0018] Step 5: Constructing a temporal frame sequence: Analyze all cell images in the image sequence, taking the cell image to be analyzed as the current frame. Using the current frame as the center, dynamically select its neighboring frames as reference frames to construct a temporal frame sequence for analysis. For the first frame, select the first three frames to construct the temporal frame; for the second frame, select the first four frames to construct the temporal frame; for the last frame, select the last three frames to construct the temporal frame; for the second-to-last frame, select the last four frames to construct the temporal frame; for other frames in the middle of the image sequence, select the current frame, the two frames before it, and the two frames after it, for a total of five frames to construct the temporal frame.

[0019] Step 6: Calculation of temporal difference features: For the current frame, perform multi-dimensional difference feature calculations on it and each reference frame in the temporal frame sequence constructed in Step 5. The multi-dimensional difference features include direct pixel difference features, gradient magnitude difference features, and high-frequency change difference features. After obtaining the multi-dimensional difference features, normalize and weightedly fuse each difference feature to obtain a fused feature value. Then, perform statistical calculations in the front-end fan-ring detection region and the back-end fan-ring detection region respectively to obtain the front-end fan-ring detection region score and the back-end fan-ring detection region score. Then, use temporal weights to weightedly fuse the scores of multiple reference frames to obtain the final comprehensive score of the front-end fan-ring detection region and the back-end fan-ring detection region of the current frame.

[0020] Step 7, Confidence Calculation and Movement Direction Determination: The final comprehensive scores of the front and back fan ring detection regions of the current frame are normalized, and the confidence is calculated based on the relative relationship between the scores of the two regions. By comparing the confidence with the set threshold, the location of the flagella is inferred. Combined with the biological propulsion mode of the flagellates, the movement direction of the flagellate cells in the current frame is finally determined.

[0021] Furthermore, the experimental apparatus described in step 1 includes a pressure source, a pressure pump, a storage bottle, a flow meter, a microchannel, and a waste bottle connected in sequence. The pressure source and the pressure pump, and the pressure pump and the storage bottle are all connected via connecting pipes. The storage bottle and the flow meter, the flow meter and the inlet of the microchannel, and the outlet of the microchannel and the waste bottle are all connected via conduits. The apparatus also includes a microscope and a high-speed camera connected in series. The high-speed camera is mounted below the stage of the microscope, and the microchannel is set on the stage of the microscope. The pressure source is used to provide the base pressure for the pressure pump; the pressure pump is used to output a set constant pressure. The system employs a pressure difference to create an internal and external pressure difference within a storage bottle. The storage bottle is sealed and contains an algal cell solution. Under the pressure difference, the algal cell solution is transported through a conduit to a microchannel, forming a shear flow. A flow meter is used to monitor the flow rate pumped in by the pressure pump in real time. The microchannel is used to form the shear flow. A waste bottle is used to collect waste algal cell solution flowing out of the microchannel outlet. The microscope uses a high-power objective lens (40-100x) to magnify microscopic images, with the standard being that algal cell flagella are visible to the naked eye. A high-speed camera is used to capture image sequences of algal cells in the shear flow.

[0022] Furthermore, the image sequence in step 1 includes at least 5 images; the effective motion trajectory refers to the fact that the flagella of the flagellated algae cells can be seen in each frame of the motion trajectory, and if flagella information is missing in 1-3 consecutive frames, it is also considered as an effective motion trajectory.

[0023] Furthermore, the Gaussian blur denoising described in step 3 specifically involves performing a convolution operation on the image using a Gaussian kernel function with a standard deviation σ = 1.0, calculated as follows:

[0024] (1)

[0025] (2)

[0026] In the formula: This represents the pixel grayscale value at coordinates (i,j) in the original image; This represents the pixel grayscale value at coordinates (x, y) after Gaussian blur processing. Represents the Gaussian kernel function;

[0027] The background estimation and removal are specifically as follows: the background is estimated using morphological opening operation, and then the background estimate is subtracted from the original image and negative values ​​are truncated to zero;

[0028] The background estimation formula is as follows:

[0029] (3)

[0030] In the formula: This represents the image after Gaussian blurring. S represents the morphological opening operation; S represents the structuring element.

[0031] The formula for background removal is as follows:

[0032] (4)

[0033] In the formula: This represents the pixel grayscale value of the original image at coordinates (x, y); This represents the estimated pixel grayscale value of the background at coordinates (x, y); This represents the pixel grayscale value at coordinates (x, y) of the image after background removal;

[0034] The contrast enhancement specifically involves: using adaptive histogram equalization to enhance local contrast, and then mapping the image intensity to the entire dynamic range through linear contrast adjustment.

[0035] The CLAHE calculation formula is as follows:

[0036] (5)

[0037] In the formula: This represents a grayscale image after background removal. This represents the enhanced image after CLAHE processing; ClipLimit represents the contrast limiting threshold parameter; NumTiles represents the image block parameters;

[0038] The formula for contrast adjustment is as follows:

[0039] (6)

[0040] In the formula: This represents the enhanced image after CLAHE processing; This represents a grayscale image after contrast adjustment; the `imadjust` function performs linearization of the image's grayscale values, transforming the input image... The grayscale value range is automatically stretched to the full [0,1] dynamic range.

[0041] Furthermore, the symmetrical region that expands to both sides at a certain angle in step 4 specifically refers to a symmetrical region that expands to both sides by 20° to 90°.

[0042] Furthermore, the calculation formula for the direct pixel difference feature mentioned in step 6 is as follows:

[0043] (7)

[0044] In the formula: This represents the grayscale value of the current frame at coordinates (x, y); This represents the grayscale value of the reference frame at coordinates (x, y); This represents the value of the direct pixel difference at coordinates (x, y);

[0045] The formula for calculating the gradient magnitude difference feature is as follows:

[0046] (8)

[0047] (9)

[0048] (10)

[0049] In the formula: , This represents the gradient components in the x and y directions of the current frame; , This represents the gradient components of the reference frame in the x and y directions; , This represents the gradient magnitude of the current frame and the reference frame at coordinates (x, y); This represents the value of the gradient magnitude difference at coordinates (x, y);

[0050] The calculation formula for the high-frequency variation difference characteristics is as follows:

[0051] (11)

[0052] (12)

[0053] (13)

[0054] In the formula: * indicates convolution operation; This represents the Gaussian kernel function with a standard deviation of σ = 2. This represents the grayscale value of the current frame at coordinates (x, y); This represents the grayscale value of the reference frame at coordinates (x, y); This represents the value of the high-frequency component of the current frame at coordinates (x, y); This represents the value of the high-frequency components of other frames at coordinates (x, y); This represents the value of the high-frequency variation difference at coordinate (x, y);

[0055] The formula for calculating the fused feature value by normalizing and weighting each differential feature is as follows:

[0056] (14)

[0057] (15)

[0058] (16)

[0059] (17)

[0060] In the formula: , , They are respectively , , The normalized value; max(D) represents the maximum value in the difference matrix D; The feature values ​​are fused; the weighting coefficients of 0.3, 0.4, and 0.3 correspond to the importance of direct pixel differences, gradient magnitude differences, and high-frequency change differences, respectively.

[0061] The formula for calculating the score of the front-end sector ring detection area is as follows:

[0062] (18)

[0063] The formula for calculating the score of the rear sector ring detection area is as follows:

[0064] (19)

[0065] In the formula: sort_descending(A) means sorting set A in descending order; N_front and N_rear represent the total number of pixels in the front fan ring detection region and the back fan ring detection region, respectively; ceil(0.25×N) means taking the smallest integer not less than 0.25×N; S_front and S_rear represent the temporal difference scores of the front fan ring detection region and the back fan ring detection region, respectively.

[0066] The calculation formula for weighted fusion of the scores of multiple reference frames using temporal weights is as follows:

[0067] (20)

[0068] (twenty one)

[0069] (twenty two)

[0070] In the formula: d_i represents the time distance between the i-th reference frame and the current frame; σ represents the standard deviation of the Gaussian weight function, which is window length / 3; w_i represents the normalized weight of the i-th reference frame; S_f_final and S_r_final represent the scores of the front-end fan-ring detection region and the back-end fan-ring detection region after weighted fusion, respectively.

[0071] Furthermore, the specific process of normalization in step 7 is as follows: Normalization is performed on the front-end sector ring detection region score S_f_final and the back-end sector ring detection region score S_r_final, calculated using the following formula:

[0072] (twenty three)

[0073] (twenty four)

[0074] (25)

[0075] In the formula: This represents the sum of the scores for the front and rear sector ring detection areas; , These represent the normalized scores for the front and rear sector detection regions, respectively.

[0076] The specific process for calculating the confidence level is as follows: The confidence level is calculated based on the normalized scores of the front and rear fan-ring detection regions. When the confidence level is lower than a set threshold, the whiplash position of the frame is determined to be unknown. The calculation formula is as follows:

[0077] (26)

[0078] In the formula: C represents the confidence level, and its range is [0,1];

[0079] The specific process for inferring the location of the flagellum is as follows: Based on the normalized scores of the front and rear fan-ring detection areas, the flagellum location is inferred; the judgment rule is as follows: if… and When the flagellum is at its tip, the flagellum is located at the front. and When the flagellum is located at the posterior end, the flagellum is positioned at the posterior end.

[0080] The specific process for determining the movement direction of the flagellated algal cells in the current frame is as follows: based on the flagella position and the direction of the cell's long axis, combined with the propulsion mode of the flagellated algae, the movement direction of the flagellated algae is determined; for algal cells that rely on the flagella to generate pulling force, the movement direction is consistent with the direction of the pulling force generated by the flagella; for algal cells that rely on the flagella to generate thrust, the movement direction is opposite to the direction of the pulling force generated by the flagella.

[0081] The beneficial effects of this invention are: the method described in this invention overcomes the technical problem of being unable to identify flagellates in dynamic water environments due to poor imaging quality and loss of instantaneous information. Through temporal dynamic feature analysis and biological motion characteristic determination, it achieves high-precision, automated, and continuous measurement of the movement direction of flagellates. Specifically, this is reflected in the following aspects:

[0082] (1) Significantly enhances the signal-to-noise ratio of weak flagellar motion features. Traditional methods rely on directly identifying the static morphology of flagella in a single frame image. However, due to the small size of the flagellar structure and its low contrast with the background, it is easily affected by noise in dynamic flow fields, leading to recognition failure. This invention analyzes the pixel changes, gradient changes, and high-frequency fluctuations of a specific detection area around the cell in a continuous time series to separate and amplify the dynamic swinging features of the flagella from the static background noise, thereby indirectly but robustly inferring the active position of the flagella and achieving a significant improvement in the signal-to-noise ratio.

[0083] (2) Effectively overcomes the problem of flagellar information loss in single or multiple frames of images. In shear flow fields, the rapid rotation of cells often causes flagella to momentarily detach from the focal plane, resulting in blurring in single or multiple consecutive frames of images, rendering traditional single-frame analysis methods ineffective. This invention utilizes the temporal continuity of the movement direction of flagellates, meaning that their posture and flagellar position do not undergo abrupt changes. When flagellar information is missing in a certain frame, this method can intelligently infer the most likely position of the flagella in the missing frame by analyzing the temporal difference signals excited by flagellar activity in adjacent frames. This analytical capability based on effective information from previous and subsequent frames ensures the continuity and data integrity of movement direction judgment under dynamic flow conditions, solving the problem of judgment interruption caused by the loss of instantaneous information in traditional methods.

[0084] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments. Attached Figure Description

[0085] Figure 1 This is a schematic diagram of the method flow described in this invention;

[0086] Figure 2 This is a schematic diagram of the experimental setup.

[0087] Figure 3 This is the original image sequence recorded in the experiment of Example 1;

[0088] Figure 4 This refers to the image sequence after cropping and alignment processing in Example 1;

[0089] Figure 5 The images shown are those before (left) and after (right) image preprocessing in Example 1.

[0090] Figure 6 This is an image showing the segmentation of the cell detection region in Example 1;

[0091] Figure 7 The timing frame constructed in Example 1;

[0092] Figure 8 This is a time-series difference heatmap from Example 1;

[0093] Figure 9 This is a schematic diagram showing the results of measuring the direction of algal cell movement in Example 1;

[0094] Figure 10 This is a statistical distribution diagram of the direction of algal cell movement measured in Example 1.

[0095] Figure 2 In the middle: 1. Pressure pump; 2. Air pressure source; 3. Liquid storage bottle; 4. Flow meter; 5. Microchannel; 6. Luer connector; 7. Tubing; 8. Waste liquid bottle. Detailed Implementation

[0096] This invention discloses a method for measuring the direction of movement of dinoflagellates under flowing conditions, such as... Figure 1 As shown, the method includes the following steps:

[0097] Step 1: Image Sequence Acquisition and Processing: An image sequence (at least 5 images) recording the effective movement trajectory of flagellated algae cells was acquired using the experimental setup. An effective movement trajectory is defined as one where the flagella of the flagellated algae are visible in every frame. Even if flagella information is missing for 1-3 consecutive frames, it is still considered a valid movement trajectory. The image sequence was then processed using the Trackmate open-source program to identify the algal cell bodies, fit them into ellipses, and save the cell ellipse parameters, including the centroid coordinates, the length of the major and minor axes, and the direction of the major axis.

[0098] like Figure 2 As shown, the experimental apparatus includes a pressure source 2, a pressure pump 1, a storage bottle 3, a flow meter 4, a microchannel 5, and a waste bottle 8 connected in sequence. The pressure source and the pressure pump, and the pressure pump and the storage bottle are connected by connecting pipes. The storage bottle and the flow meter, the flow meter and the inlet of the microchannel, and the outlet of the microchannel and the waste bottle are all connected by conduits 7. The inlet and outlet of the microchannel are connected to the conduits 6 by Luer connectors. The apparatus also includes a microscope and a high-speed camera connected in series. The high-speed camera is mounted below the stage of the microscope, and the microchannel is set on the stage of the microscope. The air pressure source provides the base pressure for the pressure pump; the pressure pump outputs a set constant pressure to create a pressure difference between the inside and outside of the storage bottle; the storage bottle is sealed and contains algal cell solution. Under the action of the internal and external pressure difference, the algal cell solution is transported through a conduit to the microchannel to form a shear flow; a flow meter is used to monitor the flow rate pumped in by the pressure pump in real time; the microchannel is the area where the algal solution flows, working with the pressure pump to create a stable and precise shear flow; a waste bottle is used to collect the algal cell waste solution flowing out of the microchannel outlet; the microscope uses a high-power objective lens (generally 40~100x objective lens, with the standard being that the algal cell flagella can be seen with the naked eye) to achieve microscopic image magnification; a high-speed camera is used to capture image sequences of algal cells in the shear flow.

[0099] The specific steps for obtaining image sequences using the experimental setup are as follows:

[0100] (1) Prepare an algal cell solution of a certain concentration, with the algal cell solution concentration being 1×10⁻⁶. 4 ~1×10 6 Connect and arrange instruments and equipment according to the experimental setup diagram between the number of samples per milliliter, and check the airtightness of the apparatus.

[0101] (2) Adjust the position of the microchannel so that the center of the camera's image frame corresponds to the center of the microchannel (the observation area is in the middle of the channel's inlet and outlet), while the upper and lower walls of the microchannel are parallel to the long border of the image frame.

[0102] (3) Determine the average flow rate according to the experimental conditions, set the pressure in the storage bottle using a pressure pump, and continuously adjust the pressure according to the real-time monitoring of the flow rate by the flow meter until the measured flow rate remains stable and consistent with the required average flow rate, so as to create a stable shear flow.

[0103] (4) After the flow stabilizes, start the high-speed camera and record and save the image sequence at a frame rate of 50 Hz. The specific frame rate is determined according to the experimental flow rate. The higher the flow rate, the higher the frame rate. Generally, the standard is that the direction of algal cell movement does not jump between adjacent images. In order to statistically analyze the distribution of algal cell movement direction, the experiment needs to be repeated many times until there are thousands of effective algal cell movement trajectories in the image sequence.

[0104] Step 2, Cell Image Cropping and Alignment: Read the cell ellipse parameters (centroid coordinates, length of major and minor axes, direction of major axis) output by Trackmate, set a standard image size (the size should be large enough to encompass the entire cell), and translate the centroid of each cell to the center of the image to achieve spatial alignment of the cell images; this eliminates analysis errors caused by the different positions of cells in the image. This step lays the foundation for subsequent time-series analysis in a unified coordinate system.

[0105] Step 3: Cell Image Preprocessing: The cropped and aligned cell images are preprocessed, primarily to suppress noise and enhance useful information, providing clear input data for subsequent time-series difference analysis. Image preprocessing includes Gaussian blur noise reduction, background estimation and removal, and contrast enhancement.

[0106] ① Gaussian Blur Denoising: A Gaussian kernel with a standard deviation σ = 1.0 is used to perform convolution operations on the image. The aim is to suppress high-frequency noise while preserving the detailed features of the flagella. The calculation formula is as follows:

[0107] (1)

[0108] (2)

[0109] In the formula: This represents the pixel grayscale value at coordinates (i,j) in the original image; This represents the pixel grayscale value at coordinates (x, y) after Gaussian blur processing. This represents the Gaussian kernel function.

[0110] ② Background estimation and removal: The background is estimated by morphological opening operation (the structuring element is a disk with a radius of 10 pixels). Then, the background estimate is subtracted from the original image and negative values ​​are truncated to zero to eliminate uneven background interference.

[0111] The background estimation formula is as follows:

[0112] (3)

[0113] In the formula: This represents the image after Gaussian blurring. represents the morphological opening operation; S represents the structuring element.

[0114] The formula for background removal is as follows:

[0115] (4)

[0116] In the formula: This represents the pixel grayscale value of the original image at coordinates (x, y); This represents the estimated pixel grayscale value of the background at coordinates (x, y); This represents the pixel grayscale value at coordinates (x, y) of the image after background removal.

[0117] ③ Contrast Enhancement: Adaptive histogram equalization is used to enhance local contrast, and then linear contrast adjustment is used to map the image intensity to the entire dynamic range to highlight weak flagellar signals.

[0118] The CLAHE calculation formula is as follows:

[0119] (5)

[0120] In the formula: This represents a grayscale image after background removal. This represents the enhanced image after CLAHE processing; ClipLimit represents the contrast limiting threshold parameter; NumTiles represents the image block parameter.

[0121] The formula for contrast adjustment is as follows:

[0122] (6)

[0123] In the formula: This represents the enhanced image after CLAHE processing; This represents a grayscale image after contrast adjustment; the `imadjust` function performs linearization of the image's grayscale values, transforming the input image... The grayscale value range is automatically stretched to the full [0,1] dynamic range.

[0124] Step 4: Create fan-ring detection regions: For the preprocessed cell image, create fan-ring detection regions centered on the cell centroid at both ends of the cell's long axis. The fan-ring detection region is defined as: an inner circle radius that avoids the cell body, an outer circle radius that covers the flagellar activity area, and a symmetrical region (covering the flagellar activity area) that extends to both sides at a certain angle (20°~90°) with the cell's long axis as the center. Thus, a front-end fan-ring detection region and a rear-end fan-ring detection region are formed at both ends of the cell's long axis.

[0125] Step 5: Construct a time-series frame: Analyze all cell images in the image sequence, take the cell image to be analyzed as the current frame, and dynamically select its neighboring frames as reference frames with the current frame as the center to construct a time-series frame sequence for analysis; for the first frame, select the first three frames to construct a time-series frame; for the second frame, select the first four frames to construct a time-series frame; for the last frame, select the last three frames to construct a time-series frame; for the second to last frame, select the last four frames to construct a time-series frame; for other frames in the middle of the image sequence, select the current frame, the two frames before it, and the two frames after it, for a total of five frames to construct a time-series frame.

[0126] Step 6: Calculation of Temporal Difference Features: For the current frame, multi-dimensional difference features are calculated between it and each reference frame in the temporal frame sequence constructed in Step 5, including direct pixel difference features, gradient amplitude difference features, and high-frequency variation difference features. This step amplifies weak flagellar motion signals by analyzing pixel changes, gradient changes, and high-frequency fluctuations in a specific detection area around the cell over a continuous period of time. Specifically, the calculation includes three dimensions of difference features: direct pixel difference features reflect direct temporal changes in image brightness, capturing obvious flagellar waving features; gradient amplitude difference features reflect temporal changes in image edges and textures, identifying subtle structural motion patterns; and high-frequency variation difference features detect rapid motion features through high-pass filtering, enhancing sensitivity to rapid flagellar waving.

[0127] After acquiring multi-dimensional differential features, each feature is normalized and weighted to obtain a fused feature value. Robust statistical calculations are then performed in both the front-end and back-end fan-loop detection regions, using the average of the top 25% of values ​​to obtain the region score, effectively suppressing noise interference and highlighting significant motion signals. Temporal weights are then used to weight and fuse the scores from multiple reference frames. These temporal weights are calculated based on a Gaussian function, with reference frames closer to the current frame receiving higher weights. This yields the final combined score for the front-end and back-end fan-loop detection regions of the current frame, providing a quantitative basis for subsequent flagellum position inference.

[0128] ① The formula for calculating direct pixel difference features is as follows:

[0129] (7)

[0130] In the formula: This represents the grayscale value of the current frame at coordinates (x, y); This represents the grayscale value of the reference frame at coordinates (x, y); This represents the value of the direct pixel difference at coordinates (x, y).

[0131] ② The formula for calculating the gradient magnitude difference characteristic is as follows:

[0132] (8)

[0133] (9)

[0134] (10)

[0135] In the formula: , This represents the gradient components in the x and y directions of the current frame; , This represents the gradient components of the reference frame in the x and y directions; , This represents the gradient magnitude of the current frame and the reference frame at coordinates (x, y); This represents the difference in gradient magnitude at coordinates (x, y).

[0136] ③ The calculation formula for high-frequency variation characteristics is as follows:

[0137] (11)

[0138] (12)

[0139] (13)

[0140] In the formula: * indicates convolution operation; This represents the Gaussian kernel function with a standard deviation of σ = 2. This represents the grayscale value of the current frame at coordinates (x, y); This represents the grayscale value of the reference frame at coordinates (x, y); This represents the value of the high-frequency component of the current frame at coordinates (x, y); This represents the value of the high-frequency components of other frames at coordinates (x, y); This represents the value of the high-frequency variation difference at coordinate (x, y).

[0141] ④ Feature weighted fusion: The three differential features are normalized and weighted and fused. The calculation formula is shown below:

[0142] (14)

[0143] (15)

[0144] (16)

[0145] (17)

[0146] In the formula: , , They are respectively , , The normalized value; max(D) represents the maximum value in the difference matrix D; The fusion feature values ​​are represented by weight coefficients of 0.3, 0.4, and 0.3, which correspond to the importance of direct pixel differences, gradient magnitude differences, and high-frequency variation differences, respectively.

[0147] ⑤ Score calculation for the sector ring detection area:

[0148] The formula for calculating the score of the front-end sector ring detection area is as follows:

[0149] (18)

[0150] The formula for calculating the score of the back-end sector ring detection area is as follows:

[0151] (19)

[0152] In the formula: sort_descending(A) means sorting set A in descending order; N_front and N_rear represent the total number of pixels in the front fan ring detection region and the back fan ring detection region, respectively; ceil(0.25×N) means taking the smallest integer not less than 0.25×N; S_front and S_rear represent the temporal difference scores of the front fan ring detection region and the back fan ring detection region, respectively.

[0153] ⑥ Temporal weighted fusion: The scores of multiple reference frames are weighted and fused, and the calculation formula is shown below:

[0154] (20)

[0155] (twenty one)

[0156] (twenty two)

[0157] In the formula: d_i represents the time distance between the i-th reference frame and the current frame; σ represents the standard deviation of the Gaussian weight function, which is window length / 3; w_i represents the normalized weight of the i-th reference frame; S_f_final and S_r_final represent the scores of the front-end fan-ring detection region and the back-end fan-ring detection region after weighted fusion, respectively.

[0158] Step 7, Confidence Calculation and Movement Direction Determination: The combined scores of the front and back fan ring detection regions of the current frame are normalized, and the confidence is calculated based on the relative relationship between the scores of the two regions. By comparing the confidence with the set threshold, the location of the flagella is inferred. Combined with the biological propulsion mode of the flagellates, the movement direction of the flagellate cells in the current frame is finally determined.

[0159] ① Score Normalization: The scores S_f_final for the front-end sector ring detection region and S_r_final for the back-end sector ring detection region are normalized. The calculation formula is as follows:

[0160] (twenty three)

[0161] (twenty four)

[0162] (25)

[0163] In the formula: This represents the sum of the scores for the front and rear sector ring detection areas; , These represent the normalized scores for the front and rear sector detection regions, respectively.

[0164] ② Confidence Calculation: The confidence score is calculated based on the normalized scores of the front and rear fan-ring detection regions. When the confidence score is lower than a set threshold (e.g., 0.1), the whiplash position of the frame is determined to be unknown. The calculation formula is as follows:

[0165] (26)

[0166] In the formula: C represents the confidence level, and its range is [0,1];

[0167] ③ Flagella position inference: Based on the normalized scores of the front and rear fan-ring detection areas, the flagella position is inferred; the judgment rules are as follows: If and When the flagellum is at its tip, the flagellum is located at the front. and When the flagellum is located at the posterior end, the flagellum is positioned at the posterior end.

[0168] ④ Determining the direction of movement: Based on the position of the flagella and the direction of the cell's long axis, combined with the propulsion mode of the flagellates, the direction of movement of the flagellates is determined; for algal cells that rely on the pulling force generated by the flagella, the direction of movement is consistent with the direction of the pulling force generated by the flagella; for algal cells that rely on the pushing force generated by the flagella, the direction of movement is opposite to the direction of the pulling force generated by the flagella.

[0169] Through extensive statistical analysis of experimental data, the distribution of the movement direction of flagellate cells can be obtained.

[0170] Example 1

[0171] This embodiment is an application example of the above method.

[0172] This embodiment uses the movement direction detection of *Heterosigma erythroplasma* (propelled by flagellar pull) in a flowing environment as an example. The microscope used in the experimental setup has a 50x objective lens, and the algal cell concentration is 1.2 × 10⁻⁶. 5 The algal cells were positioned at a density of 100 cells / mL, with an average flow rate of 100 μm / s, a microchannel width of 600 μm, and a depth of 150 μm. The movement direction of the algal cells was measured, and the distribution of movement directions was statistically analyzed based on extensive experimental data.

[0173] In this embodiment, the standard image size is set to 400×400 pixels when cropping and aligning cell images. For the preprocessed cell image, the created fan-ring detection region is a symmetrical region with an inner radius of 80 pixels and an outer radius of 120 pixels, extending 50° to both sides of the cell's long axis, with a total angle of 100° (covering the flagellar activity area).

[0174] The process and results of measuring the direction of movement of dinoflagellates in this embodiment are as follows: Figures 3-10 As shown. Among them, Figure 3 This is the original image sequence recorded in the experiment. Figure 4 This is a cropped and aligned image sequence. Figure 5 Images before (left) and after (right) preprocessing. Figure 6 The image is divided into cell detection areas (including the annular region and the front and rear fan-shaped detection areas). Figure 7 The time-series frames are constructed (the third image is the frame for the current direction of analysis, and the other images are reference frames). Figure 8 The heatmap shows the temporal differences between frame i and frames i-2, i-1, i+1, and i+2, from left to right. Figure 9 The result of measuring the direction of algal cell movement (red line indicates the direction of movement). Figure 10 A statistical distribution diagram for measuring the direction of algal cell movement.

[0175] Finally, it should be noted that the above description is only used to illustrate the technical solution of the present invention and not to limit it. Although the present invention has been described in detail with reference to the preferred arrangement, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solution of the present invention without departing from the spirit and scope of the technical solution of the present invention.

Claims

1. A method for measuring the direction of movement of dinoflagellates under flowing conditions, characterized in that, The method includes the following steps: Step 1: Image sequence acquisition and processing: Image sequences recording the effective movement trajectories of flagellate cells are acquired using the experimental setup. The image sequences are then processed to identify the cell bodies of the algae, fit them into ellipses, and save the cell ellipse parameters, including the centroid coordinates, the length of the major and minor axes, and the direction of the major axis. Step 2, Cell Image Cropping and Alignment: Read the cell ellipse parameters, set a standard image size that can encompass the entire cell, and translate the centroid of each cell to the center of the image to achieve spatial alignment of the cell image; Step 3, Cell Image Preprocessing: The cropped and aligned cell image is preprocessed, including Gaussian blur noise reduction, background estimation and removal, and contrast enhancement. Step 4: Create fan-ring detection regions: For the preprocessed cell image, fan-ring detection regions centered on the cell centroid are created at both ends of the cell's long axis. Specifically, the fan-ring detection regions are: the inner radius avoids the cell body, the outer radius covers the flagellar activity area, and symmetrical regions extend to both sides at a certain angle with the cell's long axis as the center; thus, a front-end fan-ring detection region and a rear-end fan-ring detection region are formed at both ends of the cell's long axis. Step 5: Constructing a temporal frame sequence: Analyze all cell images in the image sequence, taking the cell image to be analyzed as the current frame. Using the current frame as the center, dynamically select its neighboring frames as reference frames to construct a temporal frame sequence for analysis. For the first frame, select the first three frames to construct the temporal frame; for the second frame, select the first four frames to construct the temporal frame; for the last frame, select the last three frames to construct the temporal frame; for the second-to-last frame, select the last four frames to construct the temporal frame; for other frames in the middle of the image sequence, select the current frame, the two frames before it, and the two frames after it, for a total of five frames to construct the temporal frame. Step 6: Calculation of temporal difference features: For the current frame, perform multi-dimensional difference feature calculations on it and each reference frame in the temporal frame sequence constructed in Step 5. The multi-dimensional difference features include direct pixel difference features, gradient magnitude difference features, and high-frequency change difference features. After obtaining the multi-dimensional difference features, normalize and weightedly fuse each difference feature to obtain a fused feature value. Then, perform statistical calculations in the front-end fan-ring detection region and the back-end fan-ring detection region respectively to obtain the front-end fan-ring detection region score and the back-end fan-ring detection region score. Then, use temporal weights to weightedly fuse the scores of multiple reference frames to obtain the final comprehensive score of the front-end fan-ring detection region and the back-end fan-ring detection region of the current frame. Step 7, Confidence Calculation and Movement Direction Determination: The final comprehensive scores of the front and back fan ring detection regions of the current frame are normalized, and the confidence is calculated based on the relative relationship between the scores of the two regions. By comparing the confidence with the set threshold, the location of the flagella is inferred. Combined with the biological propulsion mode of the flagellates, the movement direction of the flagellate cells in the current frame is finally determined.

2. The method for measuring the direction of movement of flagellates under flowing conditions according to claim 1, characterized in that, The experimental apparatus described in step 1 includes a pressure source, a pressure pump, a storage bottle, a flow meter, a microchannel, and a waste bottle connected in sequence. The pressure source and the pressure pump, and the pressure pump and the storage bottle are all connected by connecting pipes. The storage bottle and the flow meter, the flow meter and the inlet of the microchannel, and the outlet of the microchannel and the waste bottle are all connected by conduits. The apparatus also includes a microscope and a high-speed camera connected in series. The high-speed camera is mounted below the stage of the microscope, and the microchannel is set on the stage of the microscope. The pressure source is used to provide a base pressure for the pressure pump. The pressure pump is used to output a set constant pressure to create an internal and external pressure difference in the storage bottle. The storage bottle is a closed system containing an algal cell solution. Under the action of the internal and external pressure difference, the algal cell solution is transported to the microchannel through the conduit to form a shear flow. The flow meter is used to monitor the flow rate pumped in by the pressure pump in real time; The microchannels are used to form a shear flow; The waste liquid bottle is used to collect the waste solution of algal cells flowing out of the microchannel outlet; The microscope employs a high-power objective lens to magnify microscopic images. The high-power objective lens is 40 to 100 times magnification, with the standard being that algal cell flagella can be seen with the naked eye. The high-speed camera is used to capture image sequences of algal cells in shear flow.

3. The method for measuring the direction of movement of flagellates under flowing conditions according to claim 1, characterized in that, The image sequence in step 1 includes at least 5 images; the effective motion trajectory refers to the presence of flagella of flagellated algae cells in each frame of the motion trajectory. If flagella information is missing in 1-3 consecutive frames, it is also considered as an effective motion trajectory.

4. The method for measuring the direction of movement of flagellates under flowing conditions according to claim 1, characterized in that, The Gaussian blur denoising step 3 specifically involves performing a convolution operation on the image using a Gaussian kernel function with a standard deviation σ = 1.

0. The calculation formula is as follows: （1） （2） In the formula: This represents the pixel grayscale value at coordinates (i,j) in the original image; This represents the pixel grayscale value at coordinates (x, y) after Gaussian blur processing. Represents the Gaussian kernel function; The background estimation and removal are specifically as follows: the background is estimated using morphological opening operation, and then the background estimate is subtracted from the original image and negative values ​​are truncated to zero; The background estimation formula is as follows: （3） In the formula: This represents the image after Gaussian blurring. S represents the morphological opening operation; S represents the structuring element. The formula for background removal is as follows: （4） In the formula: This represents the pixel grayscale value at coordinates (x, y) in the original image; This represents the estimated pixel grayscale value of the background at coordinates (x, y); This represents the pixel grayscale value at coordinates (x, y) of the image after background removal; The contrast enhancement specifically involves: using adaptive histogram equalization to enhance local contrast, and then mapping the image intensity to the entire dynamic range through linear contrast adjustment. The CLAHE calculation formula is as follows: （5） In the formula: This represents a grayscale image after background removal. This represents the enhanced image after CLAHE processing; ClipLimit represents the contrast limiting threshold parameter; NumTiles represents the image block parameters; The formula for contrast adjustment is as follows: （6） In the formula: This represents the enhanced image after CLAHE processing; This represents a grayscale image after contrast adjustment; the `imadjust` function performs linearization of the image's grayscale values, adjusting the input image... The grayscale value range is automatically stretched to the full [0,1] dynamic range.

5. The method for measuring the direction of movement of flagellates under flowing conditions according to claim 1, characterized in that, The symmetrical region that expands to both sides at a certain angle in step 4 specifically refers to the symmetrical region that expands to both sides by 20° to 90°.

6. The method for measuring the direction of movement of flagellates under flowing conditions according to claim 1, characterized in that, The calculation formula for the direct pixel difference feature mentioned in step 6 is as follows: （7） In the formula: This represents the grayscale value of the current frame at coordinates (x, y); This represents the grayscale value of the reference frame at coordinates (x, y); This represents the value of the direct pixel difference at coordinates (x, y); The formula for calculating the gradient magnitude difference feature is as follows: （8） （9） （10） In the formula: , This represents the gradient components in the x and y directions of the current frame; , This represents the gradient components of the reference frame in the x and y directions; , This represents the gradient magnitude of the current frame and the reference frame at coordinates (x, y); This represents the value of the gradient magnitude difference at coordinates (x, y); The calculation formula for the high-frequency variation difference characteristics is as follows: （11） （12） （13） In the formula: ∗ represents the convolution operation; This represents the Gaussian kernel function with a standard deviation of σ = 2. This represents the grayscale value of the current frame at coordinates (x, y); This represents the grayscale value of the reference frame at coordinates (x, y); This represents the value of the high-frequency component of the current frame at coordinates (x, y); This represents the value of the high-frequency components of the reference frame at coordinates (x, y); This represents the value of the high-frequency variation difference at coordinate (x, y); The formula for calculating the fused feature value by normalizing and weighting each differential feature is as follows: （14） （15） （16） （17） In the formula: , , They are respectively , , The normalized value; max(D) represents the maximum value in the difference matrix D; To fuse feature values; The weighting coefficients of 0.3, 0.4, and 0.3 correspond to the importance of direct pixel differences, gradient magnitude differences, and high-frequency change differences, respectively. The formula for calculating the score of the front-end sector ring detection area is as follows: （18） The formula for calculating the score of the rear sector ring detection area is as follows: （19） In the formula: sort_descending(A) means sorting set A in descending order; N_front and N_rear represent the total number of pixels in the front fan ring detection region and the back fan ring detection region, respectively; ceil(0.25×N) means taking the smallest integer not less than 0.25×N; S_front and S_rear represent the temporal difference scores of the front fan ring detection region and the back fan ring detection region, respectively. The calculation formula for weighted fusion of the scores of multiple reference frames using temporal weights is as follows: （20） （21） （22） In the formula: d_i represents the time distance between the i-th reference frame and the current frame; σ represents the standard deviation of the Gaussian weight function, which is window length / 3; w_i represents the normalized weight of the i-th reference frame; S_f_final and S_r_final represent the scores of the front-end fan-ring detection region and the back-end fan-ring detection region after weighted fusion, respectively.

7. The method for measuring the direction of movement of flagellates under flowing conditions according to claim 6, characterized in that, The normalization process described in step 7 is as follows: The scores S_f_final for the front-end sector ring detection region and S_r_final for the back-end sector ring detection region are normalized using the following formula: （23） （24） （25） In the formula: This represents the sum of the scores for the front and rear sector ring detection areas; , These represent the normalized scores for the front and rear sector detection regions, respectively. The specific process for calculating the confidence level is as follows: The confidence level is calculated based on the normalized scores of the front and rear fan-ring detection regions. When the confidence level is lower than a set threshold, the flagellum position is determined to be unknown. The calculation formula is as follows: （26） In the formula: C represents the confidence level, and its range is [0,1]; The specific process for inferring the location of the flagellum is as follows: Based on the normalized scores of the front and rear fan-ring detection areas, the flagellum location is inferred; the judgment rule is as follows: if… and When the flagellum is at its tip, the flagellum is located at the front. and When the flagellum is located at the posterior end, the flagellum is positioned at the posterior end. The specific process for determining the movement direction of the flagellated algal cells in the current frame is as follows: based on the flagella position and the direction of the cell's long axis, combined with the propulsion mode of the flagellated algae, the movement direction of the flagellated algae is determined; for algal cells that rely on the flagella to generate pulling force, the movement direction is consistent with the direction of the pulling force generated by the flagella; for algal cells that rely on the flagella to generate thrust, the movement direction is opposite to the direction of the pulling force generated by the flagella.

Images