A GPU neuron skeleton extraction method based on topology pruning backbone growth

Abstract

This invention discloses a GPU-based neuronal skeleton extraction method based on topology-simplified trunk growth, relating to the fields of neuroinformatics and 3D biomedical image processing. The method includes: threshold filtering and denoising of the input volume data, and retaining the largest connected component based on 3D connected component labeling to remove isolated noise; given the cell body center coordinates, determining the cell body region through spherical expansion and removing it from the foreground to highlight neural protrusions such as dendrites and axons; after obtaining the coarse skeleton, using GPU to perform topology simplification, and constructing a path scoring function combining path length and average normalized brightness, iteratively selecting and retaining the trunk path with the best score, enhancing trunk continuity and suppressing pseudo-branches caused by noise and breaks. This invention can achieve robust skeleton extraction in 3D data with noise interference and uneven brightness, and is suitable for applications such as neuronal morphology reconstruction and quantitative analysis.

Description

Technical Field

[0001] This invention relates to the fields of neuroinformatics and three-dimensional biomedical image processing technology, and in particular to a method for extracting the neuronal skeleton from three-dimensional neuronal microscopic imaging data. Background Technology

[0002] Neuronal morphology reconstruction is a crucial foundational task in brain science research and neural circuit analysis. Skeleton extraction from three-dimensional microscopic imaging data allows for the abstraction of voxel-level foreground structures into a central axis topology, facilitating subsequent tree structure analysis, morphological statistics, and visualization.

[0003] However, real-world imaging data often suffers from background noise, uneven brightness, breaks, and adhesions. Furthermore, the large volume and high intensity of cell bodies can easily lead to backbone shifts or the generation of numerous pseudo-branches. Moreover, existing methods are time-consuming to process large-scale neuronal data, failing to meet the ever-increasing demands for scanning precision and hindering subsequent research progress. Therefore, there is an urgent need for a rapid and effective method for extracting the neuronal backbone that can suppress noise interference, enhance backbone continuity, and provide fast processing. Summary of the Invention

[0004] To address the shortcomings of existing technologies, the present invention aims to provide a method for extracting the neuronal backbone. This method reduces noise and cell body interference on the backbone by filtering connected components and removing cell bodies, improves the continuity and stability of the backbone backbone by using a backbone preservation strategy, and finally accelerates the processing speed by using GPU optimization technology.

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

[0006] A GPU neuron skeleton extraction method based on topology-reduced backbone growth, applied to three-dimensional neuron microscopic imaging data, the method comprising:

[0007] S1, Data Acquisition: Under the condition of obtaining the coordinates of the cell body center, a spherical expansion is performed starting from the cell body center to determine the cell body region, and the cell body region is removed from the foreground voxel set to obtain the neural protrusion region;

[0008] S2, Data preprocessing: Threshold filtering and noise reduction are performed on the volume data to obtain a foreground voxel set;

[0009] Furthermore, 3D connected component labeling is performed on the foreground voxel set, preserving the largest connected component to remove isolated noise;

[0010] S3, Topology simplification: A coarse skeleton is generated based on the neural protrusion region, and the coarse skeleton is topology simplified on the GPU.

[0011] S4, Trunk scoring and growth: A path scoring function is constructed on the coarse skeleton. The path scoring function is related to at least the path length and the average normalized brightness of the voxels on the path. Then, the trunk growth process is accelerated by GPU. Based on the path scoring function, the trunk path with the best score is iteratively selected in the coarse skeleton and retained to obtain the target neuron skeleton.

[0012] Preferably, step S1 includes the following steps:

[0013] S11, Obtain cell body coordinates. This method uses manual annotation to obtain the cell body center coordinates. In the interface for 3D neuron volume data, the operator records the cell body center coordinates by specifying voxel indices / spatial coordinates. This serves as the input parameter for subsequent cell body spherical expansion to determine the cell body region and perform knockout processing.

[0014] Preferably, step S1 includes the following steps:

[0015] S12, cell body coordinate removal, obtaining cell body center coordinates After that, with Starting with the foreground voxel set, a spherical expansion is performed to determine the cell body region, and this region is then removed from the foreground to avoid interference from the cell body on subsequent skeletal trunk growth and branch determination. Let the foreground set after threshold segmentation and connected component filtering be... The candidate regions of the cell body are defined as follows:

[0016]

[0017] in To remove the radius from the cell body, when using iterative expansion, from Starting from the beginning, expand outwards layer by layer until the stopping condition is met. Then, cell body knockout is performed:

[0018]

[0019] in The set of neural protrusions after cell body removal is used for subsequent coarse skeleton generation, topology simplification and trunk growth steps. The definition of \ is as follows: A\B refers to the new set consisting of all elements that belong to set A but not to set B.

[0020] Preferably, step S2 includes the following steps:

[0021] S21, Threshold filtering: Gray-level threshold filtering is applied to the input 3D neuron microscopic imaging data to initially separate the foreground and background of the neurons. Let the gray-level function of the volume data be... ,in Represents voxel coordinates, with a preset threshold selected. (To preserve the fine protrusions of neurons as much as possible, a relatively low threshold is usually chosen.) The volume data is binarized to obtain the foreground mask. :

[0022]

[0023] in Represents the foreground voxels of candidate neurons. The background voxel is represented by the threshold-filtered foreground mask, which serves as the input for subsequent connected component filtering and noise removal.

[0024] Preferably, step S2 includes the following steps:

[0025] S22, Connected region labeling, used to obtain a binary foreground mask after threshold filtering. Then, for all those that satisfy The voxels are used to perform 3D connected component labeling to distinguish different foreground connected components and remove discrete noise. Under the preset 3D neighborhood connectivity criterion, the foreground voxel set is decomposed into several sets of non-overlapping connected components. , of which Both exist connected paths consisting of foreground voxels that allow them to reach each other; calculate the voxel count for each connected component. The connected component with the largest number of voxels is selected as the foreground component of the target neuron.

[0026]

[0027] Subsequently retained And delete the remaining connected components, thereby removing isolated small connected components and random noise;

[0028] The resulting foreground set is denoted as It is used in subsequent cell body removal and scaffold extraction steps.

[0029] Preferably, step S3 includes the following steps:

[0030] S31, Voxel Neighborhood Determination: During the topology simplification stage, for any candidate model voxel... (Foreground voxel), first determine whether it satisfies the neighborhood constraint of "deletable voxel" based on the connectivity of its 3D neighborhood. Two voxels in 3D space... Neighborhood relations are determined by Manhattan distance Defined as: when It is 6-connected. It is 18-connected. For 26-connectivity, define accordingly neighborhood set Subsequently Within the neighborhood of , the connected components of "model (foreground)" and "background" are judged respectively. Only when A voxel is marked as removable only if it meets all three of the following core conditions:

[0031] In its - Neighborhood There is one and only one inside. - Connected model components;

[0032] In its - Neighborhood There is one and only one inside. - A connected background component, and this background component is connected to... for - Connected;

[0033] At least one - Connected background neighbors, to adapt to GPU parallel judgment, further adopt... of The above neighborhood determination is achieved by template matching within the neighborhood:

[0034] In any sub-iteration in either direction, if The neighborhood matches the template corresponding to that direction (or its rotation around the boundary normal axis). The form of the template is marked as deletable; the template can be represented by a Boolean equation as follows:

[0035]

[0036] in For neighboring voxels, Its logic is NOT. These represent the sets of neighbor indices that are required to be "0", "1", and "at least one 1" in the template, respectively.

[0037] Preferably, step S3 includes the following steps:

[0038] S32, GPU-accelerated optimization: To improve the processing efficiency of the topology simplification stage on high-resolution 3D voxel data, this method parallelizes and optimizes memory access for "template matching / Boolean equation judgment of boundary voxels" on the GPU. First, a bit-encoding compression strategy is adopted:

[0039] along The shaft will be continuous The foreground / background states of voxels are packaged into a 32-bit unsigned integer, allowing a single CUDA thread to process one package unit at a time. This transforms the complex branching decisions per voxel into a single and efficient bitwise Boolean equation calculation, significantly improving computational density and reducing if-else branch overhead.

[0040] Preferably, in step S32,

[0041] Its bit encoding can be expressed as: for a fixed With the beginning ,definition:

[0042]

[0043] in Foreground indicator function, The Bit correspondence Voxel state, followed by rapid acquisition and caching of neighborhood data:

[0044] for Direction of neighbors, through Obtained by performing a displacement operation;

[0045] for Directional neighbors are obtained by reading adjacent packed integers and temporarily storing the obtained neighborhood data in the GPU's high-speed shared memory for subsequent template matching and unified calculation of Boolean equations, thereby improving memory access merging and reducing global memory access latency.

[0046] Preferably, step S4 includes the following steps:

[0047] S41, Core score calculation: After topological simplification to obtain the rough skeleton, the rough skeleton is considered as a connected structure composed of voxel points. Any two voxels are defined... and A path is formed by the sequence of connected voxels along the rough skeleton. ,in Indicates the first position on the path Individual points, in order to simultaneously characterize the topological dominance and imaging reliability of a path, for path Calculate the main score function:

[0048]

[0049] in This is the topological length of the path (i.e., the number of voxel points contained in the path). The average normalized brightness of the path:

[0050]

[0051] in, Voxel representation Normalized brightness values ​​in the original image. Brightness weighting coefficient (preferred default value) ),when The time score degenerates into a score based solely on topological length;

[0052] when At this time, the brighter path will receive a higher score, thus increasing the priority of structures with better imaging quality in the selection of the trunk while maintaining the topological prior that "the trunk is usually longer".

[0053] Preferably, step S4 includes the following steps:

[0054] S42, GPU-accelerated optimization: During the trunk growth stage, to ensure efficient execution of the longest and brightest trunk path search even with a large number of coarse skeleton voxels, the coarse skeleton is represented as a sparse graph structure with voxels as nodes (nodes are skeleton voxels, edges are neighborhood connectivity relationships). The cumulative calculation of path scores, optimal path selection for candidate endpoints, and trunk backtracking marking are parallelized on the GPU, using skeleton points near the cell body center as the root node set. Perform a hierarchical traversal of the skeleton graph and for each node Calculate the cumulative trunk score from the root to the node, and let the node brightness weight be:

[0055]

[0056] in This represents the normalized brightness of the node voxels in the original image. The brightness weight coefficient is then used to determine the distance from the root to the node. The cumulative score can be recursively derived as follows:

[0057]

[0058] in To satisfy the predecessor neighbor set of topological hierarchy constraints, and to synchronously record the predecessor pointer when the maximum value is obtained. This recursion can be executed in parallel on the GPU in units of "layers":

[0059] Each node in the same layer The calculations are independent of each other and are suitable for being completed in parallel using thread blocks;

[0060] Set of endpoints / leaf nodes The optimal backbone endpoint selection is achieved using parallel reduction computation:

[0061]

[0062] This yields the endpoint of the highest-rated main path, and then... pointer from Backtracking to the root set After marking the main path, nodes on that path are added to the "retained main path set." In the next iteration, the remaining structure continues to branch or expand the main path. To reduce GPU global memory access overhead, node adjacency information is stored in a compact, compressed sparse row adjacency list format. With cumulative value The architecture employs a structure array layout, and the maximum value search for endpoint scores combines intra-block shared memory reduction with atomic updates to achieve a high-throughput backbone selection and growth process.

[0063] The beneficial effects of this invention are as follows:

[0064] This invention reduces isolated noise and improves the robustness of skeleton extraction by filtering for the largest connected component.

[0065] This invention reduces the interference of cell bodies on centerline positioning by eliminating cell body regions, thereby reducing false branches;

[0066] This invention improves backbone continuity and suppresses erroneous connections caused by breaks and noise by combining a backbone preservation strategy with path length and brightness information.

[0067] This invention improves the speed of neuron skeleton extraction and processing through various GPU optimization techniques. Attached Figure Description

[0068] Figure 1 This is a flowchart of the overall algorithm of the present invention;

[0069] Figure 2 This is an example of the neuronal skeleton obtained by the method of the present invention. Detailed Implementation

[0070] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.

[0071] Example 1

[0072] The core of this invention lies in the skeleton extraction and backbone preservation of three-dimensional neuronal volume data. This is achieved using input volume data... For example, Figure 1 As shown, a method for extracting the neuronal skeleton based on three-dimensional neuronal microscopic imaging mainly includes the following steps:

[0073] S1, Cell body region identification and removal:

[0074] S1 is used to exclude cell body regions before scaffold extraction, preventing cell body clumps from interfering with subsequent topological simplification and trunk growth. The input is three-dimensional neuronal microscopic imaging volume data, with the volume data grayscale denoted as... ,in Let be the voxel coordinates; let the set of foreground voxels after threshold filtering and connected component filtering be . .

[0075] S11, Input cell body center coordinates:

[0076] Manual annotation is used as an input parameter to directly provide the cell body center coordinates, denoted as:

[0077]

[0078] in The voxel coordinates corresponding to the cell body center selected by the operator in the in vivo data are used as the center point parameters for subsequent cell body region removal.

[0079] S12, Cell body region identification and removal:

[0080] Using the input cell body center coordinates Centered on a spherical neighborhood, the cell body removal region is determined, and voxels within this region are removed from the foreground set. Specifically, let the cell body removal radius be . The cell body region is defined as:

[0081]

[0082] And perform cell body knockout to obtain the neural protrusion foreground set:

[0083]

[0084] in To exclude the cell body, the foreground voxel set is used as input for subsequent coarse skeleton generation, topology simplification, and trunk growth steps, so that the output skeleton can more accurately reflect the centerline structure of dendrites and axons.

[0085] S2, Threshold Filtering and Connected Component Selection:

[0086] S2 is used to preprocess the original three-dimensional volumetric data of neurons to preserve the foreground structure of neurons as much as possible and remove isolated noisy connected components, providing clean input for subsequent topology simplification and backbone growth. Since the signal-to-noise ratio of neuron images is often spatially uneven, a fixed threshold that is too high may accidentally delete weak dendrites / axons, while a threshold that is too low will leave residual noise. Therefore, this embodiment uses a manually set lower fixed threshold to prioritize the preservation of foreground points; residual noise will be further cleaned up through iterative deletion in the subsequent topology simplification stage.

[0087] S21, Threshold Filtering:

[0088] Let the grayscale of the input 3D volume data be... ,in Voxel coordinates. A fixed threshold set by the user is given. (Take the lower value) Binarize the volume data to obtain the foreground mask. :

[0089]

[0090] in Represents the foreground voxels of candidate neurons. This represents the background voxels. By setting a low threshold, we aim to preserve the slender processes of neurons after binarization.

[0091] S22, Connectivity Labeling and Filtering:

[0092] Even after threshold filtering, a small number of noise points may still exist, typically appearing as isolated small voxel connected regions, which can introduce erroneous branches during skeleton extraction. Considering that this algorithm processes 3D image data of a single neuron, this embodiment uses a binary foreground... Perform 3D connected component labeling to obtain the set of foreground connected components. .

[0093] The 3D component connectivity labeling algorithm can be summarized as follows:

[0094] On a 3D binary skeleton voxel grid, connected components are marked according to the specified neighborhood connectivity; the algorithm first establishes a label array of the same shape as the skeleton for each voxel, and unvisited positions are initially set to -1;

[0095] Then, the grid is scanned sequentially. When a foreground voxel is encountered that has not yet been labeled, a breadth-first search is started with that voxel as the seed: the seed is labeled with the current component number, and all its unlabeled neighbors that are within the boundary, are also foreground, and are not yet labeled are expanded in sequence using a queue, and they are assigned the same number.

[0096] After all the components are labeled, the number is incremented by 1, and the search continues for the next unlabeled foreground voxel until the entire volume data has been traversed.

[0097] Each foreground voxel in the final labeled array is assigned an integer number to its associated connected component. The number of voxels in each connected component is then calculated. Only retain the largest connected component:

[0098]

[0099] And delete all Small connected components are selected to eliminate isolated noisy connected components. The filtered set of foreground voxels is denoted as . This serves as input for subsequent steps.

[0100] S3, Topology Simplification Methods and GPU Optimization:

[0101] S3 is used to further refine the neuron foreground voxels obtained in step two into a topology-preserving coarse skeleton, and is accelerated through GPU parallel processing. This embodiment uses a 3D iterative deletion algorithm to complete the topology simplification, and performs bitwise-level optimization on the GPU for "template matching / Boolean equation judgment of boundary voxels".

[0102] S31, Topology Simplification:

[0103] 3D connectivity convention: Let the voxel coordinates be... Manhattan distance Define three-dimensional adjacency relationships: It is 6-connected. It is 18-connected. It is 26-connected.

[0104] Deletable voxel determination: For any model voxel Deletion is only permitted if its removal does not alter the image's topology; this "removable voxel" must simultaneously satisfy the three core conditions C3.1–C3.3 given in the paper:

[0105] C3.1: In its - Neighborhood There is one and only one inside. - Connected model components;

[0106] C3.2: In its - Neighborhood There is one and only one inside. - A connected background component, and this background component is connected to... for - Connected;

[0107] C3.3: At least one - Connected background neighbors.

[0108] Six-sub-iteration boundary order strategy: To ensure the safety of parallel deletion, a complete iteration is divided into six sub-iterations in six directions. Only boundary voxels facing the corresponding direction (e.g., up / down / east / west / south / north) are processed. Each sub-iteration deletes only the set of voxels that meet the conditions.

[0109] 3×3×3 template matching and rotation equivalence: In each sub-iteration, by checking voxels of Does the neighborhood match any template in that direction (and the template rotated about the boundary normal axis)? (The morphology after deletion) is used to mark the removable voxels.

[0110] Boolean equation representation of a template: For ease of computation, a template can be written in Boolean equation form. For a given template... Its determination can be summarized as follows:

[0111]

[0112] in For neighboring voxels, Logical NOT; These are the neighbor index sets that require voxels to be "0", "1", and "at least one 1" in the template, respectively.

[0113] Termination condition: The algorithm peels away removable voxels layer by layer from the model surface; when there are no voxels left to be removed in a complete iteration (all six sub-iterations have been executed), the iteration terminates, and the set of remaining voxels is the topologically simplified centerline rough skeleton.

[0114] S32, GPU acceleration optimization:

[0115] To efficiently execute the above "template Boolean judgment + iterative deletion" on the GPU, this embodiment transforms the batch judgment of voxel states into bitwise operations and optimizes neighborhood access:

[0116] Bit-code compression: Encoding voxel data along... Axis compression allows a 32-bit integer to be stored contiguously for 32 seconds. The state of the directional voxel (0 for background, 1 for model); allocate a thread on the GPU to process a 32-bit integer, thereby processing 32 voxels in parallel at once, transforming the deletion decision from complex if-else logic into a single, highly optimizable Boolean equation, and implementing it with bitwise operations such as AND / OR / NOT / XOR.

[0117] Neighborhood data access optimization: Template determination requires accessing 26 neighborhoods. For the packaged 32 voxels, Directional neighbors are quickly obtained by performing a bitwise shift operation on a 32-bit integer. and Directional neighbors are obtained by reading adjacent 32-bit integers. To improve memory access efficiency, edge-to-edge access is used. Axis encoding is used to ensure that threads within the same warp access contiguous global memory addresses; and the acquired neighborhood data is temporarily stored in GPU shared memory for reuse in subsequent Boolean equation calculations.

[0118] S4, Stemline Score and Growth:

[0119] In step S4, a hierarchical backbone preservation strategy of "from coarse to fine, with priority given to the backbone" is adopted on the coarse skeleton (connected structure composed of voxel points) obtained in step S3: first, the longest global path is determined by two breadth-first search (BFS) as the backbone, and then the secondary branches with the best quality are iteratively discovered and retained based on the determined backbone; at the same time, the multiple BFS graph traversal processes are accelerated in parallel on the GPU in the form of a double-buffered front queue to reduce traversal overhead and improve overall efficiency.

[0120] Step S4 includes the following:

[0121] Path definition and backbone scoring, for any two voxels on the rough skeleton and The voxel sequence obtained by connecting the inter-skeleton backbone is defined as a path:

[0122]

[0123] in For the first on the path Individual element.

[0124] The topological length of a path is defined as:

[0125]

[0126] The average normalized brightness of the path is defined as:

[0127]

[0128] in It represents voxels Normalized brightness values ​​in the original image.

[0129] Calculate the path scoring function based on topology length and brightness information:

[0130]

[0131] in , This is a brightness weighting coefficient (the default value can be 0.5 in the examples), used to increase the priority of paths with better imaging quality during the preservation process based on the topological prior that "the trunk is usually longer".

[0132] Branch initialization:

[0133] During the initialization phase, two BFS operations are performed starting from any voxel point in the rough skeleton: the first BFS is used to find the farthest voxel point as one end of the path; the second BFS is performed starting from this farthest point to determine the farthest point at the other end; and the longest path is obtained by backtracking using the predecessor information recorded during the second BFS. The longest path is used as the neuron backbone and added to the retention set:

[0134]

[0135] This initialization process ensures that the most important skeletal structure is preserved in a priority and complete manner, providing a stable "main trunk base" for subsequent branch iterations.

[0136] Iterative branching growth:

[0137] In the main trunk Once determined, the iteration phase begins. Round iteration with the current hold set As a "growth substrate", from All voxel points are treated as the same frontier set and BFS is performed to search for skeleton voxel points that have not yet been preserved and form a set of candidate branch paths.

[0138] To avoid retaining excessively short and insignificant branches, a length threshold constraint is imposed on candidate paths during the search process: when a candidate path is found to be shorter than the threshold for this round... When this happens, the expansion in that direction is stopped, thereby preventing short branches from entering the candidate set.

[0139] Subsequently, the candidate paths that satisfy the length constraint are calculated one by one. And select only the single path with the highest score. Add to the reserved set:

[0140]

[0141] This ensures that the "most valuable" secondary branches are retained in each round, avoiding redundancy caused by blindly retaining all branches.

[0142] Meanwhile, this embodiment employs an adaptive threshold strategy, causing the threshold to automatically decrease with each branch level: The length threshold for a round is defined as one-third of the length of the longest candidate path found in that round.

[0143]

[0144] in For the first The maximum length of the candidate path is determined by this strategy. First-level branches closer to the trunk have a higher threshold because they are longer, while higher-level branches farther from the trunk have a lower threshold because they are shorter, thus achieving a more balanced retention ratio among branches at different levels.

[0145] The iterative process terminates when any of the following conditions are met: the maximum length of a candidate path in a given round is less than the global minimum length threshold. or retain the set It covers the entire set of skeleton voxels, or reaches the preset maximum number of iterations. To prevent excessive iteration.

[0146] GPU parallel breadth-first search (BFS) optimization:

[0147] To accelerate the numerous BFS operations in the above steps, this embodiment employs a double-buffered front queue architecture on the GPU to implement layer-synchronous parallel BFS: two integer arrays (each with the capacity of the total number of voxels) are pre-allocated in the device memory as front queue buffers, serving as the "current front" and the "next front," respectively; a parallel expansion kernel function is launched in each iteration, dynamically configuring the thread grid based on the number of current front nodes, with each thread responsible for expanding the neighbors of one front node. This design strictly separates reads and writes: all threads read only from the current front, while newly discovered nodes are uniformly written to the next front, thereby eliminating read-write dependencies.

[0148] When writing to the next leading edge, the write position is allocated through an atomic counter on the device: after a thread successfully marks a new node as visited, it atomically increments the counter, and the return value serves as the unique index of the node in the queue, enabling large-scale thread conflict-free concurrent writing; after expansion, the next leading edge is changed to the new current leading edge through pointer swapping, without copying the queue contents; the final value of the counter serves as the number of nodes in the next layer, used to determine the termination and configuration of the next round of thread size.

[0149] During the branch discovery phase, to avoid affecting the retained set Executing a single-source BFS at each starting point would lead to excessive kernel function startups and redundant accesses. This embodiment employs a multi-source parallel BFS: All retained nodes are copied in batches to the initial frontier queue, and their distances are uniformly initialized to 0. This is equivalent to starting from multiple source points simultaneously to expand the wavefront, thereby merging multiple independent searches into a single parallel operation, significantly reducing the number of kernel function startups and redundant calculations.

[0150] In addition, to prevent nodes that have been retained in the branch search from being repeatedly selected as candidate path endpoints, a parallel invalidation kernel function is started before each round of branch discovery iteration to traverse voxels and mark the distance values ​​of retained nodes as invalid values ​​so that they are automatically excluded in the subsequent selection of farthest points.

[0151] Example 2

[0152] To demonstrate the effectiveness of this invention compared with other neuron skeleton extraction methods, i.e., the experimental results of different neurons and models, precision, recall, and F-score were selected as evaluation indicators, and the quantitative results are shown in Table 1.

[0153] Table 1

[0154]

[0155] The quantitative experimental results (in milliseconds) of runtime for these datasets are shown in Table 2. This invention outperforms other classic models on the complex difficulty models in Table 2, achieving the shortest runtime. It is also very fast on other difficulty models (easy or medium difficulty). The Fastmarching, Advanced, and MOST methods all fail to run successfully. Figure 2 As shown, Figure 2 This is an example of the neuronal skeleton obtained by the method of the present invention.

[0156] Table 2

[0157]

[0158] In addition to its advantages of better performance and faster speed, this invention also performs very stably on different models, thus having wide applicability in practical use.

[0159] The embodiments described above are merely examples of several implementation methods of this application, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the invention patent. It should be noted that those skilled in the art can make various modifications and improvements without departing from the concept of this application, and these all fall within the protection scope of this application.

Claims

1. A method for extracting the backbone of a GPU neuron based on topologically simplified backbone growth, characterized in that, The method, applied to three-dimensional neuronal microscopic imaging data, includes: S1, Data Acquisition: Under the condition of obtaining the coordinates of the cell body center, a spherical expansion is performed starting from the cell body center to determine the cell body region, and the cell body region is removed from the foreground voxel set to obtain the neural protrusion region; S2, Data preprocessing: Threshold filtering and denoising are performed on the volume data to obtain the foreground voxel set, and three-dimensional connected component labeling is performed on the foreground voxel set to retain the largest connected component to remove isolated noise; S3, Topology simplification: A coarse skeleton is generated based on the neural protrusion region, and the coarse skeleton is topology simplified on the GPU. S4, Trunk scoring and growth: A path scoring function is constructed on the coarse skeleton. The path scoring function is related to at least the path length and the average normalized brightness of the voxels on the path. Then, the trunk growth process is accelerated by GPU. Based on the path scoring function, the trunk path with the best score is iteratively selected in the coarse skeleton and retained to obtain the target neuron skeleton.

2. The method for extracting the GPU neuron backbone based on topologically simplified backbone growth according to claim 1, characterized in that, Step S1 includes the following steps: S11, Obtain cell body coordinates. This method uses manual annotation to obtain the cell body center coordinates. In the interface for 3D neuron volume data, the operator records the cell body center coordinates by specifying voxel indices / spatial coordinates. This serves as the input parameter for subsequent cell body spherical expansion to determine the cell body region and perform knockout processing.

3. The method for extracting the GPU neuron backbone based on topologically simplified backbone growth according to claim 1, characterized in that, Step S1 includes the following steps: S12, cell body coordinate removal, obtaining cell body center coordinates After that, with Starting with the foreground voxel set, a spherical expansion is performed to determine the cell body region, and this region is removed from the foreground. Let the foreground set after threshold segmentation and connected component filtering be... The candidate regions of the cell body are defined as follows: in To remove the radius from the cell body, when using iterative expansion, from Starting from the beginning, expand outwards layer by layer until the stopping condition is met. Then, cell body knockout is performed: in This is a foreground set of neural protrusions after cell body removal, used for subsequent rough skeleton generation, topology simplification, and trunk growth steps.

4. The method for extracting the GPU neuron backbone based on topologically simplified backbone growth according to claim 1, characterized in that, Step S2 includes the following steps: S21, Threshold filtering: Gray-level threshold filtering is applied to the input 3D neuron microscopic imaging data to initially separate the foreground and background of the neurons. Let the gray-level function of the volume data be... ,in Represents voxel coordinates, with a preset threshold selected. Binarize the volume data to obtain the foreground mask : in Represents the foreground voxels of candidate neurons. The background voxel is represented by the threshold-filtered foreground mask, which serves as the input for subsequent connected component filtering and noise removal.

5. The method for extracting the GPU neuron backbone based on topologically simplified backbone growth according to claim 1, characterized in that, Step S2 includes the following steps: S22, Connected region labeling, used to obtain a binary foreground mask after threshold filtering. Then, for all those that satisfy The voxels are labeled with 3D connected components. Under a preset 3D neighborhood connectivity criterion, the foreground voxel set is decomposed into several sets of non-overlapping connected components. , of which Both exist connected paths consisting of foreground voxels that make them reachable. Calculate the voxel count for each connected component. The connected component with the largest number of voxels is selected as the foreground component of the target neuron. Subsequently retained And delete the remaining connected components, thereby removing isolated small connected components and random noise; the resulting foreground set is denoted as It is used in subsequent cell body removal and scaffold extraction steps.

6. The method for extracting the GPU neuron backbone based on topologically simplified backbone growth according to claim 1, characterized in that, Step S3 includes the following steps: S31, Voxel Neighborhood Determination: During the topology simplification stage, for any candidate model voxel... First, based on the connectivity within its three-dimensional neighborhood, it is determined whether it satisfies the neighborhood constraint of a "deletable voxel". Two voxels in three-dimensional space... Neighborhood relations are determined by Manhattan distance Defined as: when It is 6-connected. It is 18-connected. For 26-connectivity, define accordingly neighborhood set Subsequently Within the neighborhood of "model" and "background", the connected components of "model" and "background" are judged respectively. Only when... A voxel is marked as removable only if it meets all three of the following core conditions: In its - Neighborhood There is one and only one inside. - Connected model components; In its - Neighborhood There is one and only one inside. - A connected background component, and this background component is connected to... for - Connected; At least one - Connected background neighbors, to adapt to GPU parallel judgment, further adopt... of The above neighborhood determination is achieved by template matching within the neighborhood: In any sub-iteration in either direction, if If the neighborhood matches the template corresponding to that direction, it is marked as deletable; The template can be represented by a Boolean equation as follows: in For neighboring voxels, Its logic is NOT. These represent the sets of neighbor indices in the template that are required to be "0", "1" and "at least one 1", respectively.

7. The method for extracting the GPU neuron backbone based on topologically simplified backbone growth according to claim 1, characterized in that, Step S3 includes the following steps: S32 performs parallelization and memory access optimization on the GPU for "template matching / Boolean equation judgment of boundary voxels". First, a bit-encoding compression strategy is adopted: along... The shaft will be continuous The foreground / background states of voxels are packaged into a 32-bit unsigned integer, allowing a single CUDA thread to process one package unit at a time, transforming the complex branching judgments per voxel into a single and efficient bitwise Boolean equation calculation.

8. The method for extracting the GPU neuron backbone based on topologically simplified backbone growth according to claim 1, characterized in that, In step S32 Its bit encoding can be expressed as: for a fixed With the beginning ,definition: in Foreground indicator function, The Bit correspondence Voxel state, followed by rapid acquisition and caching of neighborhood data: for Direction of neighbors, through Obtained by performing a displacement operation; for The directional neighbors are obtained by reading the adjacent packed integers, and the obtained neighborhood data is temporarily stored in the GPU's high-speed shared memory for subsequent template matching and unified calculation of Boolean equations.

9. The method for extracting the GPU neuron backbone based on topologically simplified backbone growth according to claim 1, characterized in that, Step S4 includes the following steps: S41, Core score calculation: After topological simplification to obtain the rough skeleton, the rough skeleton is considered as a connected structure composed of voxel points. Any two voxels are defined... and A path is formed by the sequence of connected voxels along the rough skeleton. ,in Indicates the first position on the path Individual points, in order to simultaneously characterize the topological dominance and imaging reliability of a path, for path Calculate the main score function: in The topological length of the path. The average normalized brightness of the path: in, Voxel representation Normalized brightness values ​​in the original image. For brightness weighting coefficients, when The scoring degenerates into a scoring based solely on topological length; when At that time, the brighter path will receive a higher score.

10. The method for extracting the GPU neuron backbone based on topologically simplified backbone growth according to claim 1, characterized in that, Step S4 includes the following steps: S42, GPU-accelerated optimization: During the trunk growth stage, the coarse skeleton is represented as a sparse graph structure with voxels as nodes. The cumulative calculation of path scores, the selection of optimal paths for candidate endpoints, and the trunk backtracking marking are parallelized on the GPU. The skeleton points near the cell body center are used as the root node set. Perform a hierarchical traversal of the skeleton graph and for each node Calculate the cumulative trunk score from the root to the node, and let the node brightness weight be: in This represents the normalized brightness of the node voxels in the original image. The brightness weight coefficient is then used to determine the distance from the root to the node. The cumulative score can be recursively derived as follows: in To satisfy the predecessor neighbor set of topological hierarchy constraints, and to synchronously record the predecessor pointer when the maximum value is obtained. This recursion can be executed in parallel on the GPU in units of "layers": Each node in the same layer The calculations are independent of each other and are suitable for being completed in parallel using thread blocks; Set of endpoints / leaf nodes The optimal backbone endpoint selection is achieved using parallel reduction computation: This yields the endpoint of the highest-rated main path, and then... pointer from Backtracking to the root set After marking the main path, nodes on that main path are added to the "retained main path set". In the next iteration, the remaining structure continues to branch or expand the main path. Node adjacency information is stored in a compact, compressed sparse row adjacency list format. Node brightness... With cumulative value The algorithm uses a structure array layout, and the maximum value search for endpoint scores combines intra-block shared memory reduction with atomic updates.

Images