Photoacoustic blood vessel segmentation method and device based on three-way maximum intensity projection intersection
Patent Information
- Application Number
- CN202611441437.8
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-09-16
- Publication Date
- 2026-10-27
AI Technical Summary
然而,上述分割方法存在明显的局限性:首先,三维手动标注的效率极低,逐层手动圈定区域耗时可达数十分钟至数小时,且操作者间差异性大,导致结果的重复性差;其次,自动分割算法高度依赖大量的标注数据,且由于主要针对CT(Computed Tomography,计算机断层扫描)/MRI(Magnetic Resonance Imaging,磁共振成像)数据进行训练,光声图像特有的声学伪影与非均匀信噪比会使其迁移精度显著下降;再者,现有分割方法均未涉及三向MIP掩膜在三维空间中的逻辑交集操作,无法直接用于无训练数据场景下的确定性几何分割
[0021]本申请提供的基于三向最大强度投影交集的光声血管分割方法和装置,沿三个正交轴向对原始三维光声图像进行最大强度投影生成二维图像,巧妙地将复杂的立体空间能量压缩聚集于高信噪比的二维平面,不仅从源头上规避了直接进行三维多维运算带来的巨大计算复杂度和内存消耗,随后在二维投影上确定目标区域并生成二维二值掩膜,将极易受空间伪影干扰的三维空间分割降维转化为边界更清晰、技术更成熟的二维图像处理,有效避开了直接三维分割的空间多义性。通过将二维掩膜反向拉伸回三维空间并执行三维空间逻辑交集处理来构建三维血管分割掩膜,在立体空间中引入了独特的正交几何形态学约束机制,如同拉起了一张多维视线一致性的高精度拦截网,以极低的计算成本瞬间剔除那些仅在单一视角表现为血管而实质是离焦伪影、声学反射的非血管结构,从而快速、精准地截取输出高保真的初步分割图像,实现了空间拓扑结构的几何粗筛。在此基础上,通过两阶段背景纯化对初步分割图像进行深度定制化精筛,第一阶段手段紧扣真实血管在空间中必然连续延伸、绝不孤立悬浮的解剖学先验形态特征,通过三维局部连通特征分析,将交集阶段逃逸但空间孤立的微小噪声体素强制置零,达成拓扑结构层面的去毛刺;第二阶段手段则紧扣光声成像光散与声散射的物理信号强度分布特征,以最大响应值的预设比例实施全局背景硬截断,消除血管外围雾霾状的弥散底噪。最终将带有严重伪影和散射底噪的原始光声数据,转化为具备纯黑零背景、边缘锐利、拓扑高度连通的最终三维血管图像。
Smart Images

Figure CN122989719A_ABST
Abstract
Description
Technical Field
[0001] This application relates to the field of photoacoustic vessel segmentation technology, and in particular to a photoacoustic vessel segmentation method and apparatus based on the intersection of three-dimensional maximum intensity projections. Background Technology
[0002] Photoacoustic imaging (PAI), a novel biomedical imaging technology that combines the high contrast of optics with the deep penetration of ultrasound, can visualize vascular networks at high resolution and high signal-to-noise ratio within millimeter to centimeter depths without the need for contrast agents. Taking liver vessels as a typical example, three-dimensional photoacoustic imaging can non-invasively acquire the spatial distribution and morphological characteristics of blood vessels. This has significant clinical value for applications such as tumor angiogenesis assessment, liver fibrosis grading, and surgical planning. Therefore, precise segmentation and post-processing of three-dimensional photoacoustic vessels are crucial.
[0003] Currently, in practical applications, blood vessel segmentation is mainly performed through the following methods: First, the operator manually delineates the region of interest (ROI) of the blood vessel layer by layer in three-dimensional voxel space; second, an automatic segmentation algorithm based on deep learning is used; and third, sparse annotations on maximum intensity projection (MIP) maps are used to supervise three-dimensional centerline segmentation, or two-dimensional MIP annotations are used as neural network supervision signals, and segmentation is performed offline through machine learning model training. However, the above segmentation methods have obvious limitations: First, the efficiency of three-dimensional manual annotation is extremely low, with manual delineation of the region layer by layer taking tens of minutes to several hours, and the large differences between operators lead to poor repeatability of the results; second, automatic segmentation algorithms are highly dependent on a large amount of labeled data, and since they are mainly trained on CT (Computed Tomography) / MRI (Magnetic Resonance Imaging) data, the acoustic artifacts and non-uniform signal-to-noise ratio unique to photoacoustic images will significantly reduce their transfer accuracy; third, existing segmentation methods do not involve the logical intersection operation of three-dimensional MIP masks in three-dimensional space, and cannot be directly used for deterministic geometric segmentation in scenarios without training data.
[0004] Therefore, there is an urgent need for a method to achieve efficient, accurate, and repeatable deterministic geometric segmentation of photoacoustic three-dimensional vascular data. Summary of the Invention
[0005] In view of this, this application provides a photoacoustic blood vessel segmentation method and apparatus based on the intersection of three-dimensional maximum intensity projections, so as to achieve efficient, accurate and repeatable deterministic geometric segmentation of photoacoustic three-dimensional blood vessel data.
[0006] Specifically, this application is implemented through the following technical solution:
[0007] The first aspect of this application provides a photoacoustic vessel segmentation method based on the intersection of three-dimensional maximum intensity projections, the method comprising:
[0008] A three-dimensional photoacoustic image to be processed is acquired, and the three-dimensional photoacoustic image is projected with maximum intensity along three mutually orthogonal axes to generate three corresponding two-dimensional projection images.
[0009] The target blood vessel regions are determined on the two-dimensional projection image, and three corresponding two-dimensional binary masks are generated.
[0010] The two-dimensional binary mask is extended to the three-dimensional voxel space along the corresponding orthogonal axes to construct three three-dimensional extended masks;
[0011] Perform three-dimensional spatial logical intersection processing on the three-dimensional extended mask to construct a three-dimensional blood vessel segmentation mask;
[0012] Based on the three-dimensional blood vessel segmentation mask, the blood vessel structure is extracted from the three-dimensional photoacoustic image, and a preliminary segmentation image is output;
[0013] The preliminary segmented image is subjected to two-stage background purification to output a purified three-dimensional blood vessel segmentation image; wherein, the two-stage background purification includes isolated noise removal processing based on spatial local connectivity features and background truncation processing based on global signal intensity.
[0014] A second aspect of this application provides a photoacoustic vessel segmentation device based on the intersection of three-dimensional maximum intensity projections, the device comprising a generation module, a construction module, and an output module;
[0015] The generation module is used to acquire the three-dimensional photoacoustic image to be processed, and to project the three-dimensional photoacoustic image with maximum intensity along three mutually orthogonal axes to generate three corresponding two-dimensional projection images.
[0016] The generation module is also used to determine the target blood vessel region on the two-dimensional projection image and generate three corresponding two-dimensional binary masks.
[0017] The construction module is used to extend the two-dimensional binary mask to the three-dimensional voxel space along the corresponding orthogonal axes to construct three three-dimensional extended masks.
[0018] The construction module is also used to perform three-dimensional spatial logical intersection processing on the three-dimensional extended mask to construct a three-dimensional blood vessel segmentation mask;
[0019] The output module is used to extract the vascular structure from the three-dimensional photoacoustic image based on the three-dimensional vascular segmentation mask and output a preliminary segmentation image;
[0020] The output module is also used to perform two-stage background purification on the preliminary segmented image and output a purified three-dimensional blood vessel segmentation image; wherein, the two-stage background purification includes isolated noise removal processing based on spatial local connectivity features and background truncation processing based on global signal intensity.
[0021] The photoacoustic vessel segmentation method and apparatus based on the intersection of three-dimensional maximum intensity projections provided in this application generate a two-dimensional image by projecting the original three-dimensional photoacoustic image with maximum intensity along three orthogonal axes. This cleverly compresses and concentrates the complex three-dimensional spatial energy into a two-dimensional plane with a high signal-to-noise ratio. This not only avoids the huge computational complexity and memory consumption caused by direct three-dimensional multidimensional operations from the source, but also determines the target region on the two-dimensional projection and generates a two-dimensional binary mask. This transforms the three-dimensional spatial segmentation, which is highly susceptible to spatial artifacts, into two-dimensional image processing with clearer boundaries and more mature technology, effectively avoiding the spatial ambiguity of direct three-dimensional segmentation. By stretching the two-dimensional mask back into three-dimensional space and performing three-dimensional spatial logical intersection processing to construct a three-dimensional vessel segmentation mask, a unique orthogonal geometric morphological constraint mechanism is introduced into the three-dimensional space. This is like setting up a high-precision interception net with multi-dimensional line-of-sight consistency, instantly eliminating non-vascular structures that only appear as vessels from a single viewpoint but are actually defocus artifacts or acoustic reflections with extremely low computational cost. This allows for the rapid and accurate extraction and output of a high-fidelity preliminary segmentation image, achieving geometric coarse screening of spatial topology. Building upon this foundation, a two-stage background purification process is employed to deeply customize and refine the initially segmented images. The first stage focuses on the anatomical priori morphological characteristics of real blood vessels, which are inherently continuous and never isolated in space. Through 3D local connectivity analysis, tiny noise voxels that escaped from the intersection phase but remained spatially isolated are forcibly zeroed out, achieving despiking at the topological level. The second stage focuses on the physical signal intensity distribution characteristics of photoacoustic imaging's light and sound scattering. A global hard background truncation is implemented at a preset ratio based on the maximum response value, eliminating the hazy, diffuse background noise around the blood vessels. Finally, the original photoacoustic data, burdened with severe artifacts and scattering noise, is transformed into a final 3D blood vessel image with a pure black background, sharp edges, and high topological connectivity. Attached Figure Description
[0022] Figure 1 A flowchart of the photoacoustic vessel segmentation method based on the intersection of three-dimensional maximum intensity projections provided in Embodiment 1 of this application;
[0023] Figure 2 This is a schematic diagram of the photoacoustic blood vessel segmentation device based on the intersection of three-dimensional maximum intensity projections provided in Embodiment 2 of this application. Detailed Implementation
[0024] Exemplary embodiments will now be described in detail, examples of which are illustrated in the accompanying drawings. When the following description relates to the drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with this application.
[0025] The terminology used in this application is for the purpose of describing particular embodiments only and is not intended to be limiting of the application. The singular forms “a,” “the,” and “the” used herein are also intended to include the plural forms unless the context clearly indicates otherwise. It should also be understood that the term “and / or” as used herein refers to and includes any and all possible combinations of one or more of the associated listed items.
[0026] It should be understood that although the terms first, second, third, etc., may be used in this application to describe various information, such information should not be limited to these terms. These terms are only used to distinguish information of the same type from one another. For example, without departing from the scope of this application, first information may also be referred to as second information, and similarly, second information may also be referred to as first information. Depending on the context, the word "if" as used herein may be interpreted as "when," "when," or "in response to determination."
[0027] The following specific embodiments are given to illustrate the technical solution of this application in detail.
[0028] Figure 1 This is a flowchart of the photoacoustic vessel segmentation method based on the intersection of three-directional maximum intensity projections provided in Embodiment 1 of this application. Please refer to... Figure 1 The method provided in this embodiment may include:
[0029] S101. Obtain the three-dimensional photoacoustic image to be processed, and project the three-dimensional photoacoustic image with maximum intensity along three mutually orthogonal axes to generate three corresponding two-dimensional projection images.
[0030] Specifically, a 3D photoacoustic image refers to a three-dimensional digital matrix (data volume) reconstructed using the "photoacoustic effect," reflecting the spatial distribution of light-absorbing substances within biological tissues. In a computer, a 3D photoacoustic image is represented as a three-dimensional matrix composed of "voxels (3D pixels)." Each voxel has coordinates and a corresponding grayscale response value (signal intensity). In this image, blood vessels, rich in hemoglobin, exhibit extremely strong photoacoustic signals, appearing as "high-brightness structures," while surrounding ordinary tissues have weaker signals, appearing as a "dark background."
[0031] Maximum Intensity Projection (MIP) is a computer graphics algorithm used for visualization and dimensionality reduction of 3D volumetric data. Imagine a beam of parallel light penetrating a 3D blood vessel cube, with a sheet of white paper (the projection surface) placed on the other side of the cube. As the light beam passes through a row of voxels along a selected axis (e.g., the Z-axis), the algorithm searches for the grayscale values of all voxels along this path and extracts only the value of the brightest (highest intensity) voxel, recording it at the corresponding position on the white paper. A 2D projected image refers to a planar second-order pixel matrix (i.e., a regular 2D image) produced after processing a 3D photoacoustic image using MIP. This is essentially three "flattened" 2D photographs taken from three completely perpendicular viewpoints (typically corresponding to the X, Y, and Z axes of a Cartesian coordinate system).
[0032] In specific implementation, the three-dimensional photoacoustic image is projected with maximum intensity along the first axis to generate a first two-dimensional projection image located in the plane formed by the second and third axes; the three-dimensional photoacoustic image is projected with maximum intensity along the second axis to generate a second two-dimensional projection image located in the plane formed by the first and third axes; the three-dimensional photoacoustic image is projected with maximum intensity along the third axis to generate a third two-dimensional projection image located in the plane formed by the first and second axes; the first, second, and third axes are orthogonal to each other; the maximum intensity values of the first, second, and third two-dimensional projection images are obtained respectively, and the upper limit of the display dynamic range of each two-dimensional projection image is set to the product of the maximum intensity value of the corresponding two-dimensional projection image and a preset contrast coefficient.
[0033] Specifically, the maximum pixel intensity value is retrieved along each projection path of the entire 3D photoacoustic image along the first axis (X-axis) and mapped onto a 2D plane (YZ plane) formed by the second and third axes to generate a first 2D projection image. The maximum pixel intensity value is then retrieved along each projection path of the entire 3D photoacoustic image along the second axis (Y-axis) and mapped onto a 2D plane (XZ plane) formed by the first and third axes to generate a second 2D projection image. The maximum pixel intensity value is then retrieved along each projection path of the entire 3D photoacoustic image along the third axis (Z-axis) and mapped onto a 2D plane (XY plane) formed by the first and second axes to generate a third 2D projection image. All planar pixels in the first 2D projection image are traversed, and the maximum intensity value is determined and recorded. The same process is repeated for the second and third 2D projection images. The maximum intensity value of the first two-dimensional projected image is multiplied by a preset contrast coefficient (e.g., 0.6), and the calculated product is set as the upper limit of the display dynamic range of the first two-dimensional projected image. The maximum intensity value of the second two-dimensional projected image is multiplied by the preset contrast coefficient, and the calculated product is set as the upper limit of the display dynamic range of the second two-dimensional projected image. The maximum intensity value of the third two-dimensional projected image is multiplied by the preset contrast coefficient, and the calculated product is set as the upper limit of the display dynamic range of the third two-dimensional projected image.
[0034] S102. Determine the target blood vessel regions on the two-dimensional projection image and generate three corresponding two-dimensional binary masks.
[0035] Specifically, the target vascular region refers to the region of interest containing the main vascular network structure, delineated by the operator in a two-dimensional projection image. Although maximum intensity projection (MIP) has flattened and extracted the blood vessels from the three-dimensional data, some irrelevant tissue background, acoustic artifacts, or other non-target structures that are not within the scope of this application may still remain in the two-dimensional projection image.
[0036] A two-dimensional binary mask is a special type of planar digital image (or pixel matrix). The size of a two-dimensional binary mask is exactly the same as a two-dimensional projected image, but its pixels no longer have complex grayscale levels (such as rich brightness variations from 0 to 255). Each pixel has only two possible values: either 0 or 1. A value of 1 (active / white) corresponds to pixels inside the target blood vessel region, representing blood vessels that need to be preserved. A value of 0 (inactive / black) corresponds to all background pixels outside the target blood vessel region, representing impurities that need to be filtered out.
[0037] In specific implementation, the system receives a region selection command input by the user, collects the target edge trajectory or vertex coordinate sequence on each of the two-dimensional projection images, constructs a closed boundary covering the target blood vessel structure on each of the two-dimensional projection images, and determines the target blood vessel region; it then determines the geometric inside-outside position relationship between each pixel in each of the two-dimensional projection images and the corresponding closed boundary, activates the pixels determined to be inside the closed boundary, and keeps the remaining pixels inactive, generating three two-dimensional binary masks.
[0038] Specifically, the system receives region selection commands (such as polygon clicks, including but not limited to rectangular selection, lasso tools, and brush annotations) input by the user for each 2D projected image. In response to the user's region selection command, the system collects the mouse movement edge trajectory drawn by the user in real time on each 2D projected image, or records a series of polygon vertex coordinate sequences generated by continuous user clicks. The collected edge trajectories or vertex coordinate sequences are connected end-to-end in the current plane to construct a geometrically closed boundary enclosing the target blood vessel structure on each 2D projected image, thus defining the target blood vessel region. For each 2D projected image, a boundary-inside / outside determination algorithm (such as the inpolygon criterion) is invoked to calculate and determine the spatial geometric position relationship between the coordinates of each planar pixel in the image and the geometrically closed boundary (i.e., whether the point is inside or outside the boundary). Based on the determination result, the pixel values are logically assigned: pixels determined to be inside the geometrically closed boundary are activated (their pixel values are uniformly set to an active state value, such as 1); all other pixels determined to be outside the geometrically closed boundary remain inactive (their pixel values are uniformly set to a background value, such as 0). After determining and binarizing all pixels of the three two-dimensional projection images, three independent two-dimensional binary masks are generated and output. For example, the three generated two-dimensional binary masks are: M1(y, z) (YZ projection), M2(x, z) (XZ projection), and M3(x, y) (XY projection).
[0039] S103. Extend the two-dimensional binary mask to the three-dimensional voxel space along the corresponding orthogonal axes to construct three three-dimensional extended masks.
[0040] Specifically, three-dimensional voxel space refers to a three-dimensional digital matrix composed of countless regularly arranged tiny cubic units stacked in three dimensions of space. Ordinary digital images (two-dimensional images) are planes composed of pixels arranged in a horizontal and vertical grid, while voxels are the extension of pixels in three-dimensional space, i.e., three-dimensional pixels or small cubes. Each voxel in three-dimensional voxel space has a unique and definite three-dimensional spatial coordinate (x, y, z), and internally stores the corresponding photoacoustic signal intensity value.
[0041] A 3D extended mask refers to a stereoscopic binary mask constructed in 3D voxel space by stretching a 2D binary mask along an axis perpendicular to it. In the preceding steps, the 3D image is flattened into a 2D projection along a certain axis, while the expansion is the reverse of this process. For example, suppose the first 2D binary mask is a blood vessel contour map located in the XY plane. This planar contour is copied along the perpendicular Z-axis throughout the entire 3D voxel space. In this way, in the generated 3D space, regardless of how the Z-axis coordinate changes, the shape of any XY section is exactly the same as the original 2D mask.
[0042] In specific implementation, the two-dimensional binary mask generated by projection onto the first orthogonal axis is copied and expanded along the first orthogonal axis until the expanded data dimension is consistent with the number of voxels of the three-dimensional photoacoustic image along the first orthogonal axis, thus constructing a first three-dimensional expanded mask; the two-dimensional binary mask generated by projection onto the second orthogonal axis is copied and expanded along the second orthogonal axis until the expanded data dimension is consistent with the number of voxels of the three-dimensional photoacoustic image along the second orthogonal axis, thus constructing a second three-dimensional expanded mask; the two-dimensional binary mask generated by projection onto the third orthogonal axis is copied and expanded along the third orthogonal axis until the expanded data dimension is consistent with the number of voxels of the three-dimensional photoacoustic image along the third orthogonal axis, thus constructing a third three-dimensional expanded mask.
[0043] Specifically, a planar two-dimensional binary mask previously generated corresponding to the projection along the first orthogonal axis (X-axis) is obtained. Using this two-dimensional binary mask as a reference plane, data copying and stretching are performed along the first orthogonal axis. Copying continues until the dimension of the expanded stereo data is exactly the same as the number of voxels in the original three-dimensional photoacoustic image along the first orthogonal axis, at which point the first three-dimensional expanded mask is finally constructed and output. A planar two-dimensional binary mask previously generated corresponding to the projection along the second orthogonal axis is obtained. Using this two-dimensional binary mask as a reference plane, data copying and stretching are performed along the second orthogonal axis. Copying continues until the dimension of the expanded stereo data is exactly the same as the number of voxels in the original three-dimensional photoacoustic image along the second orthogonal axis, at which point the second three-dimensional expanded mask is finally constructed and output. A planar two-dimensional binary mask previously generated corresponding to the projection along the third orthogonal axis is obtained. Using this two-dimensional binary mask as a reference plane, data copying and stretching are performed along the third orthogonal axis. The process continues until the expanded 3D data dimension is exactly the same as the number of voxels in the third orthogonal axis of the original 3D photoacoustic image, and then stops. Finally, the third 3D expanded mask is constructed and output.
[0044] For example, three three-dimensional extended masks can be represented as:
[0045] Ω1(x,y,z) = repmat(M1(y,z),[X,1,1]);
[0046] Ω2(x,y,z) = repmat(M2(x,z),[1,Y,1]);
[0047] Ω3(x,y,z) = repmat(M3(x,y),[1,1,Z]);
[0048] Wherein, Ω1 is the first three-dimensional extended mask; Ω2 is the second three-dimensional extended mask; Ω3 is the third three-dimensional extended mask; M1, M2, and M3 are three two-dimensional binary masks.
[0049] S104. Perform three-dimensional spatial logical intersection processing on the three-dimensional extended mask to construct a three-dimensional blood vessel segmentation mask.
[0050] Specifically, the 3D blood vessel segmentation mask is a 3D binary data volume (voxel matrix) that is ultimately constructed by performing logical intersection (intersection by voxel) operations on three orthogonal stereo masks in 3D space, and accurately locks the 3D geometric boundary of the target blood vessel.
[0051] In practice, the previously constructed first, second, and third 3D extended masks are obtained, and it is confirmed that the 3D voxel dimensions (length, width, and height) of the three masks are completely consistent. A 3D blank matrix (all zeros) with the same size as the 3D photoacoustic image is created and initialized in memory to store the finally constructed 3D blood vessel segmentation mask. A 3D spatial coordinate loop is started, traversing the voxel points with the same 3D spatial coordinates (x, y, z) in the first, second, and third 3D extended masks. A logical AND operation is performed on the three voxel values corresponding to the current spatial coordinates (x, y, z) (or the three corresponding voxel values are directly multiplied point by point). If the voxel values of the current spatial coordinates are all active (e.g., all values are 1) in all three 3D extended masks, then the voxel point of the corresponding coordinate in the 3D blank matrix is assigned an active value (e.g., 1). If the voxel value of the current spatial coordinates is inactive in any 3D extended mask (e.g., there is a value of 0), then the voxel point of the corresponding coordinate in the 3D blank matrix remains inactive (e.g., 0). This process continues until every coordinate point in the 3D spatial voxel array has been traversed, completing the full-image assignment of the 3D blank matrix, and finally constructing and outputting the 3D blood vessel segmentation mask.
[0052] For example, a three-dimensional blood vessel segmentation mask can be represented as:
[0053] Ωseg(x,y,z) = Ω1AND Ω2AND Ω3;
[0054] It can also be expressed as:
[0055] (Ω1 + Ω2 + Ω3) = 3;
[0056] Wherein, Ωseg is the three-dimensional blood vessel segmentation mask; Ω1 is the first three-dimensional extended mask; Ω2 is the second three-dimensional extended mask; and Ω3 is the third three-dimensional extended mask.
[0057] S105. Based on the three-dimensional blood vessel segmentation mask, extract the blood vessel structure from the three-dimensional photoacoustic image and output a preliminary segmentation image.
[0058] Specifically, the vascular structure refers to the real three-dimensional vascular network hidden in the original 3D photoacoustic image, presented by hemoglobin-rich blood targets, and its corresponding photoacoustic feature signals (grayscale amplitude data). The vascular structure includes not only the geometric topology of the blood vessels in three-dimensional space (branching, direction, diameter, spatial coordinates), but also the most original photoacoustic reconstructed grayscale values (absolute magnitude of absorber intensity) at these spatial locations. The 3D vascular segmentation mask obtained in the preceding steps is essentially just a set of purely geometric "0 and 1 switches" (it only knows where the blood vessels are, but lacks brightness information). The vascular structure to be extracted in this step is to map the regions with values of 1 in the mask back to the original 3D photoacoustic image, completely extracting the original photoacoustic signal energy (intensity variations) that truly have medical diagnostic value. The preliminary segmentation image refers to the 3D data volume output after multiplying the original 3D photoacoustic image and the 3D vascular segmentation mask voxel-by-voxel matrix (or masking), retaining only the original signal of the target blood vessel.
[0059] In specific implementation, a three-dimensional Gaussian filter is used to smooth the three-dimensional blood vessel segmentation mask, softening the boundary of the three-dimensional blood vessel segmentation mask and preserving the gray-level gradient information of the edge transition region to generate a smoothed three-dimensional blood vessel segmentation mask; the smoothed three-dimensional blood vessel segmentation mask is multiplied voxel by voxel with the three-dimensional photoacoustic image to output a preliminary segmentation image.
[0060] Specifically, a 3D Gaussian filter with a pre-defined spatial kernel size and standard deviation is invoked or constructed within the system. The input 3D blood vessel segmentation mask is fed into this filter for 3D spatial convolution calculation, and smoothing filtering is performed on the mask image. Through the blurring and diffusion effect of Gaussian filtering, the steep, hard boundaries of "0 and 1" in the original binary mask are transformed into stepped transition regions containing continuous gray-level gradient information, thereby softening the boundaries and outputting a smoothed 3D blood vessel segmentation mask. The generated smoothed 3D blood vessel segmentation mask is then fully aligned with the original 3D photoacoustic image in terms of spatial coordinates and dimensions. For each voxel coordinate point in 3D space, the voxel value of the smoothed 3D blood vessel segmentation mask at the corresponding coordinate is multiplied point-by-point with the original photoacoustic signal intensity value (gray-level value) of the 3D photoacoustic image at the same coordinate. After completing the multiplication operation for all corresponding voxel points in the entire image, the resulting new 3D volume data is used as the preliminary segmentation image.
[0061] Optionally, after outputting the preliminary segmentation image, the method further includes: performing linear normalization on the preliminary segmentation image; using a three-dimensional median filter with a set kernel size to perform three-dimensional filtering and denoising on the normalized image, and outputting a denoised three-dimensional image; using a blood vessel enhancement filter based on the Hessian matrix to perform enhancement filtering on the denoised three-dimensional image, enhancing the morphological response value of tubular blood vessel structures in space and suppressing non-tubular noise structures by analyzing the feature value structure of local second-order intensity changes in the image, and outputting a three-dimensional enhanced image; and replacing the preliminary segmentation image with the three-dimensional enhanced image to perform two-stage background purification.
[0062] Specifically, the initial segmented image is read, and its overall pixel / voxel grayscale values are linearly normalized to output a normalized image. A pre-defined 3D median filter (e.g., 7×7×7 voxels) is used to perform spatial 3D filtering on the normalized image. The median of voxel intensity within the filter kernel window is used to replace the central voxel value, thereby filtering out image noise and outputting a denoised 3D image. The denoised 3D image is then input into a Hessian matrix-based vessel enhancement filter. A Hessian matrix is constructed by calculating the second-order partial derivatives for the local region of each voxel, and this matrix is mathematically decomposed to extract the feature value structure reflecting the local second-order intensity changes of the image. Based on the calculated feature value structure, geometric morphology is determined. For regions determined to be tubular vessel structures in space, their morphological response values are proportionally enlarged or increased; for regions determined to be non-tubular noise structures, amplitude suppression and reduction are performed. After filtering, a 3D enhanced image is output. This 3D enhanced image is used as the new standard input to directly replace the initial segmented image in the original process. The replaced 3D enhanced image data stream is passed and input into the subsequent execution module for a two-stage background purification operation.
[0063] For example, in one embodiment, the three-dimensional filtering and denoising process can be expressed as:
[0064] img_denoised = medfilt3(img_norm, [7, 7, 7]);
[0065] Where img_denoised is the denoised 3D image; img_norm is the preliminary segmented image.
[0066] Enhanced filtering can be expressed as:
[0067] V = fibermetric(img_denoised, scale=8, ObjectPolarity='bright');
[0068] Where V is the 3D enhanced image; img_denoised is the denoised 3D image; scale is the scale parameter used to match the thickness (diameter) of the target blood vessel. When the size of scale matches the typical radius or diameter range of photoacoustic blood vessels in the 3D image, the morphological response value generated by the filter reaches its maximum, thereby accurately locking the blood vessel at that thickness level; ObjectPolarity is the target polarity parameter, specifying the brightness attribute (foreground color polarity) of the target structure relative to the background. Setting it to 'bright' means that the algorithm specifically identifies and enhances "bright tubular objects in a black background".
[0069] S106. Perform two-stage background purification on the preliminary segmented image and output the purified three-dimensional blood vessel segmentation image.
[0070] The two-stage background purification includes isolated noise removal based on spatial local connectivity features and background truncation based on global signal strength.
[0071] Specifically, although the above steps have filtered out large chunks of non-target tissue and visible impurities, the initial segmented image still retains two types of fatal noise in its microscopic data due to limitations of the photoacoustic reconstruction algorithm itself, acoustic scattering, and thermal noise from the system hardware. One type is spatially isolated noise outside the blood vessels, i.e., some randomly distributed, high-amplitude photoacoustic noise or local structural artifacts. The other type is continuous global diffuse background noise, i.e., low-amplitude background signals diffused throughout the three-dimensional space like "smog." If such an image with "spicules" and "smog" is directly output, it will produce significant errors when measuring blood vessel diameter, branch angle, or blood vessel density. Therefore, a two-stage background purification process is needed to address these two types of noise separately.
[0072] In specific implementation, an adaptive threshold is used to binarize the preliminary segmented image to obtain a binary image; isolated connected regions with a volume smaller than a set volume threshold are identified and removed from the binary image, and the response values corresponding to the isolated connected regions in the preliminary segmented image are forcibly set to zero, outputting the first-stage processed image; using a preset proportion of the maximum response value in the first-stage processed image as a global truncation threshold, the response values of all voxels in the first-stage processed image below the global truncation threshold are forcibly set to zero, obtaining a purified three-dimensional blood vessel segmentation image.
[0073] Specifically, the Otsu algorithm is used to calculate the adaptive threshold of the initial segmentation image. 0.5 times this adaptive threshold is taken as the final binarization threshold, and the initial segmentation image is binarized to obtain the corresponding 3D binary image. In the 3D binary image, 3D connected component labeling and searching are performed to identify each independent connected component and calculate the number of voxels contained within each connected component. The number of voxels in each connected component is compared with a preset volume threshold (2000 voxels), locking isolated connected components with a volume strictly less than 2000 voxels. The voxel coordinates of these locked isolated connected components in 3D space are indexed, and the response values (grayscale values) at the same coordinates in the initial segmentation image are forcibly set to zero. After clearing all isolated components with less than 2000 voxels in the entire image, the first-stage processed image is output. The response values of all voxels in the first-stage processed image are statistically analyzed, and the maximum response value is retrieved and recorded. This maximum response value is multiplied by a preset ratio (15%) to calculate the specific global truncation threshold (i.e., 15% of the maximum response value). The process iterates through each voxel in the first-stage processed image, comparing the current voxel's response value with the global truncation threshold one by one. If the current voxel's response value is lower than the global truncation threshold, its response value is forced to zero; if it is higher than or equal to the threshold, the original response value remains unchanged. After completing the threshold comparison and truncation processing for all voxels in the image, the purified 3D blood vessel segmentation image is finally output.
[0074] Optionally, after outputting the purified three-dimensional blood vessel segmentation image, the method further includes: performing maximum intensity projection calculations on the purified three-dimensional blood vessel segmentation image along three mutually orthogonal axes to generate corresponding original two-dimensional projection images from three perspectives; performing linear normalization processing on each of the original two-dimensional projection images; rotating each normalized two-dimensional projection image by 90 degrees to adjust and align the display coordinate system of the image; and exporting the two-dimensional projection images from each perspective after correcting the display coordinate system according to preset resolution parameters.
[0075] Specifically, the purified 3D blood vessel segmentation image is acquired. Maximum signal intensity is retrieved and extracted from voxels along each projection line of the 3D data along three mutually orthogonal axes (X-axis, Y-axis, Z-axis). The extracted maximum intensity values are compressed and mapped onto a 2D plane, generating three original 2D projection images from different viewpoints. Each generated original 2D projection image undergoes linear grayscale mapping to normalize its pixel value range to a unified display interval. Each normalized 2D projection image is then rotated 90 degrees counterclockwise or clockwise within its respective 2D plane. Through rotation transformation, the physical pixel matrix of each image is adjusted and reconstructed, completing the alignment and correction of the image display coordinate system. Using the system's preset resolution parameters (e.g., 600 DPI), the 2D projection images from each viewpoint after coordinate system correction are spatially resampled or resized to conform to the preset resolution specifications. The 2D projection images from each viewpoint after resolution adjustment are converted to a specified standard image format for data encoding and written to disk, finally outputting the corresponding 2D projection image file.
[0076] Optionally, after outputting the purified 3D blood vessel segmentation image, the method further includes: performing binarization processing on the purified 3D blood vessel segmentation image to obtain a 3D binary blood vessel volume; applying a 3D topology-preserving skeletonization algorithm to extract the blood vessel centerline skeleton; and calculating multi-parameter blood vessel morphology and topological feature indices based on the blood vessel centerline skeleton and the 3D binary blood vessel volume; the multi-parameter blood vessel morphology and topological feature indices include total blood vessel volume, total blood vessel length, average blood vessel diameter, and number of blood vessel branch points.
[0077] Specifically, a purified 3D blood vessel segmentation image is obtained. The voxel response values of this 3D image are binarized and thresholded, assigning voxels belonging to the blood vessel region a value of 1 (activated state) and voxels in the background region a value of 0 (inactive state), constructing and outputting a 3D binary blood vessel volume. Using the 3D binary blood vessel volume as input, a 3D topology-preserving skeletonization algorithm is applied. While strictly maintaining the original spatial topology and connectivity of the blood vessel, the outer boundary voxels of the 3D blood vessel are peeled and thinned layer by layer inward until the blood vessel structure is simplified to a 3D cue network of single voxel widths, finally extracting and outputting the blood vessel centerline skeleton. The extracted blood vessel centerline skeleton is spatially correlated with the 3D binary blood vessel volume, and the following four feature indicators are quantitatively calculated sequentially. Total blood vessel volume calculation: The total number of activated voxels (value 1) in the 3D binary blood vessel volume is counted and multiplied by the actual physical volume parameter corresponding to a single voxel to calculate the total blood vessel volume. Total vessel length calculation: Traverse and count the 3D spatial distances between all adjacent voxels in the vessel centerline skeleton, and sum the distances of all skeleton line segments in the entire image to calculate the total vessel length. Average vessel diameter calculation: Based on the geometric model, divide the calculated total vessel volume by the total vessel length, and convert using a scaling factor (or directly use a 3D distance transformation algorithm to calculate the average distance from the skeleton point to the nearest vessel boundary) to calculate the average vessel diameter. Number of vessel branch points calculation: Perform a neighborhood topological connectivity scan on each pixel node in the vessel centerline skeleton, identify and count the total number of intersection and bifurcation points (Nodes) connecting three or more vessel segments on the skeleton line, and thus determine the number of vessel branch points.
[0078] For example, in one embodiment, the total volume of the blood vessels can be expressed as:
[0079] Vol = sum(BW(:)) × vox_size 3 ;
[0080] Where Vol is the total volume of the blood vessel; BW is the three-dimensional binary volume of the blood vessel; sum(BW(:)) is the total number of blood vessel voxels; and vox_size is the physical side length of the voxels.
[0081] The total length of blood vessels can be expressed as:
[0082] Length = sum(skel(:)) × vox_size;
[0083] Where Length is the total length of the blood vessel; skel is the binary voxel of the blood vessel centerline skeleton; sum(skel(:)) is the total number of skeleton voxels; and vox_size is the physical side length of the voxels.
[0084] The average blood vessel diameter can be expressed as:
[0085] Diameter = 2 × mean(D(skel)) × vox_size;
[0086] Where Diameter is the average vessel diameter; skel is the binary form of the vessel centerline skeleton; D is the distance transformation matrix, representing the local radius of the vessel (voxel unit); mean is used for mean calculation; and vox_size is the physical side length of the voxel.
[0087] The number of vascular branch points can be expressed as:
[0088] BranchCount = sum(conv3(skel,ones(3,3,3)));
[0089] Where BranchCount is the number of vascular branch points; skel is the binary volume of the vascular centerline skeleton, skel > 3; conv3 is the three-dimensional spatial convolution.
[0090] The method provided in this embodiment, in its first aspect, involves stretching three two-dimensional binary masks in opposite directions along corresponding orthogonal axes into a three-dimensional voxel space, and performing a three-dimensional spatial logical intersection to construct a unified three-dimensional blood vessel segmentation mask, thereby capturing and outputting a preliminary segmented image. A multi-view spatial geometric morphological constraint mechanism is introduced, which is like setting up a precise orthogonal interception net in three-dimensional space. Utilizing the consistency of the three-dimensional line of sight, it instantly geometrically removes large-scale non-vascular structures that appear as blood vessels in a single or dual-view projection but are actually defocus artifacts or acoustic reflection artifacts. Without the need for complex deformation models or energy-intensive three-dimensional deep learning, it can quickly and accurately lock the three-dimensional geometric envelope boundary of the real three-dimensional blood vessel network with extremely low computational cost, completing a faithful restoration from two-dimensional features to three-dimensional topology.
[0091] Secondly, the first stage of local removal closely adheres to the prior morphological characteristics of real blood vessels in anatomy, which are necessarily continuous and interconnected in space, extending into a network. Through connectivity analysis, it rejects and forcibly sets to zero those tiny artifact patches that escaped during the geometric intersection stage but are isolated and suspended in space, achieving despiking at the topological level. The second stage of global truncation closely adheres to the physical distribution characteristics of photoacoustic signal intensity. Using a preset proportion of the global maximum intensity as a hard threshold, it eliminates the weak, diffuse background signal, resembling haze, formed around the blood vessels due to tissue light scattering and acoustic scattering. The two stages of purification complement each other, purifying the initial image with microscopic spikes and diffuse background into a final three-dimensional blood vessel image with absolutely pure black background, clear and sharp edges, and high topological connectivity. This provides a high signal-to-noise ratio and 100% deterministic standard data source for downstream medical clinical quantitative analysis of blood vessel morphology (such as diameter measurement and branch counting).
[0092] Thirdly, after obtaining high-purity three-dimensional vascular images, multi-parameter quantitative calculations are performed, achieving a leap from qualitative vision to precise quantitative diagnosis. The width and centerline of individual voxels are extracted using binarization and three-dimensional topology-preserving skeletonization algorithms, and then multi-parameter features are calculated by combining distance transformation. First, the activated voxels of the binary volume are statistically analyzed, and the total volume and length of the blood vessels are calculated by combining them with physical side lengths. This provides absolutely accurate physical quantitative values reflecting the filling and extension scale of microvessels, eliminating the errors of traditional subjective judgment. Second, the mean of distance transformation on the skeleton points is extracted to convert the average vessel diameter, enabling the capture of diameter distortion of each branch axis at extremely high resolution, providing highly sensitive deformation warnings for early lesions of vascular stenosis or aneurysms. Third, the number of branch points is statistically counted by scanning the skeleton nodes through three-dimensional spatial convolution, quantitatively characterizing the complexity of the vascular network from a topological perspective, providing a gold standard for the clinical grading of tumor neovascularization. In summary, combining skeleton extraction with geometric feature quantization perfectly unlocks the benefits of high signal-to-noise ratio data after upstream denoising and purification, constructing a multi-dimensional, highly deterministic objective digital diagnostic matrix, which greatly improves the scientific nature of auxiliary decision-making in clinical vascular lesion screening and efficacy evaluation.
[0093] Corresponding to the aforementioned embodiment of the photoacoustic vessel segmentation method based on the intersection of three-dimensional maximum intensity projections, this application also provides an embodiment of the photoacoustic vessel segmentation device based on the intersection of three-dimensional maximum intensity projections.
[0094] Figure 2 This is a schematic diagram of the photoacoustic vessel segmentation device based on the intersection of three-dimensional maximum intensity projections provided in Embodiment 2 of this application. Please refer to... Figure 2 The apparatus provided in this embodiment includes a generation module 210, a construction module 220, and an output module 230;
[0095] The generation module 210 is used to acquire the three-dimensional photoacoustic image to be processed, and to project the three-dimensional photoacoustic image with maximum intensity along three mutually orthogonal axes to generate three corresponding two-dimensional projection images.
[0096] The generation module 210 is also used to determine the target blood vessel region on the two-dimensional projection image and generate three corresponding two-dimensional binary masks.
[0097] The construction module 220 is used to extend the two-dimensional binary mask to the three-dimensional voxel space along the corresponding orthogonal axes to construct three three-dimensional extended masks.
[0098] The construction module 220 is also used to perform three-dimensional spatial logical intersection processing on the three-dimensional extended mask to construct a three-dimensional blood vessel segmentation mask;
[0099] The output module 230 is used to extract the vascular structure from the three-dimensional photoacoustic image based on the three-dimensional vascular segmentation mask and output a preliminary segmentation image;
[0100] The output module 230 is also used to perform two-stage background purification on the preliminary segmented image and output a purified three-dimensional blood vessel segmentation image; wherein, the two-stage background purification includes isolated noise removal processing based on spatial local connectivity features and background truncation processing based on global signal intensity.
[0101] The apparatus of this embodiment can be used to perform... Figure 1 The steps of the method embodiment shown are similar in principle and process, and will not be repeated here.
[0102] The specific implementation process of the functions and roles of each unit in the above device can be found in the implementation process of the corresponding steps in the above method, and will not be repeated here.
[0103] For the device embodiments, since they basically correspond to the method embodiments, the relevant parts can be referred to in the description of the method embodiments. The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this application according to actual needs. Those skilled in the art can understand and implement this without creative effort.
[0104] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the scope of protection of this application.
Claims
1. A photoacoustic vessel segmentation method based on the intersection of three-dimensional maximum intensity projections, characterized in that, The method includes: A three-dimensional photoacoustic image to be processed is acquired, and the three-dimensional photoacoustic image is projected with maximum intensity along three mutually orthogonal axes to generate three corresponding two-dimensional projection images. The target blood vessel regions are determined on the two-dimensional projection image, and three corresponding two-dimensional binary masks are generated. The two-dimensional binary mask is extended to the three-dimensional voxel space along the corresponding orthogonal axes to construct three three-dimensional extended masks; Perform three-dimensional spatial logical intersection processing on the three-dimensional extended mask to construct a three-dimensional blood vessel segmentation mask; Based on the three-dimensional blood vessel segmentation mask, the blood vessel structure is extracted from the three-dimensional photoacoustic image, and a preliminary segmentation image is output; The preliminary segmented image is subjected to two-stage background purification to output a purified three-dimensional blood vessel segmentation image; wherein, the two-stage background purification includes isolated noise removal processing based on spatial local connectivity features and background truncation processing based on global signal intensity.
2. The method according to claim 1, characterized in that, The preliminary segmented image undergoes two-stage background purification to output a purified 3D blood vessel segmentation image, including: The preliminary segmented image is binarized using an adaptive threshold to obtain a binary image; In the binary image, isolated connected regions with a volume smaller than a set volume threshold are identified and removed. The response value of the isolated connected regions in the preliminary segmentation image is forced to zero, and the first-stage processed image is output. Using a preset proportion of the maximum response value in the first-stage processed image as a global truncation threshold, the response values of all voxels in the first-stage processed image that are below the global truncation threshold are forced to zero, thereby obtaining a purified three-dimensional blood vessel segmentation image.
3. The method according to claim 1, characterized in that, The two-dimensional binary mask is extended into the three-dimensional voxel space along the corresponding orthogonal axes to construct three three-dimensional extended masks, including: The two-dimensional binary mask generated by projection corresponding to the first orthogonal axis is copied and expanded along the first orthogonal axis until the expanded data dimension size is consistent with the number of voxels of the three-dimensional photoacoustic image on the first orthogonal axis, thus constructing the first three-dimensional expanded mask. The two-dimensional binary mask generated by projection corresponding to the second orthogonal axis is copied and expanded along the second orthogonal axis until the expanded data dimension size is consistent with the number of voxels of the three-dimensional photoacoustic image on the second orthogonal axis, thus constructing a second three-dimensional expanded mask. The two-dimensional binary mask generated by projecting the image along the third orthogonal axis is copied and expanded along the third orthogonal axis until the expanded data dimension is consistent with the number of voxels of the three-dimensional photoacoustic image along the third orthogonal axis, thus constructing a third three-dimensional expanded mask.
4. The method according to claim 1, characterized in that, A three-dimensional photoacoustic image to be processed is acquired, and the image is projected with maximum intensity along three mutually orthogonal axes to generate three corresponding two-dimensional projection images, including: The three-dimensional photoacoustic image is projected with maximum intensity along the first axis to generate a first two-dimensional projection image located on the plane formed by the second and third axes; The three-dimensional photoacoustic image is projected with maximum intensity along the second axis to generate a second two-dimensional projection image located in the plane formed by the first axis and the third axis; The three-dimensional photoacoustic image is projected with maximum intensity along the third axis to generate a third two-dimensional projection image located in the plane formed by the first axis and the second axis; the first axis, the second axis and the third axis are orthogonal to each other. The maximum intensity values of the first two-dimensional projection image, the second two-dimensional projection image, and the third two-dimensional projection image are obtained respectively. The upper limit of the display dynamic range of each two-dimensional projection image is set to the product of the maximum intensity value of the corresponding two-dimensional projection image and the preset contrast coefficient.
5. The method according to claim 1, characterized in that, The target blood vessel regions are determined on the two-dimensional projection image, and three corresponding two-dimensional binary masks are generated, including: The system receives a region selection command input by the user, acquires the target edge trajectory or vertex coordinate sequence on each of the two-dimensional projection images, constructs a closed boundary covering the target blood vessel structure on each of the two-dimensional projection images, and determines the target blood vessel region. The geometric positional relationship between each pixel in each of the two-dimensional projection images and the corresponding closed boundary is determined one by one. Pixels that are determined to be inside the closed boundary are activated, while the remaining pixels remain inactive, generating three two-dimensional binary masks.
6. The method according to claim 1, characterized in that, Based on the aforementioned 3D blood vessel segmentation mask, the blood vessel structure is extracted from the 3D photoacoustic image, and a preliminary segmentation image is output, including: A three-dimensional Gaussian filter is used to smooth the three-dimensional blood vessel segmentation mask, softening the boundary of the three-dimensional blood vessel segmentation mask and preserving the gray-level gradient information of the edge transition region, thereby generating a smoothed three-dimensional blood vessel segmentation mask. The smoothed 3D blood vessel segmentation mask is multiplied voxel-by-voxel by the 3D photoacoustic image to output a preliminary segmentation image.
7. The method according to claim 1, characterized in that, After outputting the preliminary segmented image, the method further includes: The preliminary segmented image is linearly normalized, and a three-dimensional median filter with a set kernel size is used to perform three-dimensional filtering and denoising on the normalized image to output the denoised three-dimensional image. The denoised 3D image is enhanced by a blood vessel enhancement filter based on the Hessian matrix. By analyzing the feature value structure of local second-order intensity changes in the image, the morphological response value of tubular blood vessel structures in space is enhanced and non-tubular noise structures are suppressed, and a 3D enhanced image is output. The initial segmentation image is replaced with the three-dimensional enhanced image to perform two-stage background purification.
8. The method according to claim 1, characterized in that, After outputting the purified 3D blood vessel segmentation image, the method further includes: The purified three-dimensional blood vessel segmentation image is binarized to obtain a three-dimensional binary blood vessel volume; A three-dimensional topology-preserving skeletonization algorithm is applied to extract the vessel centerline skeleton. Based on the vessel centerline skeleton and the three-dimensional binary vessel volume, multi-parameter vessel morphology and topological feature indicators are calculated. The multi-parameter vessel morphology and topological feature indicators include the total vessel volume, total vessel length, average vessel diameter, and number of vessel branch points.
9. The method according to claim 1, characterized in that, After outputting the purified 3D blood vessel segmentation image, the method further includes: Maximum intensity projection calculations are performed on the purified three-dimensional blood vessel segmentation image along three mutually orthogonal axes to generate the original two-dimensional projection images from the corresponding three perspectives. Each of the original two-dimensional projection images is linearly normalized, and each normalized two-dimensional projection image is rotated by 90 degrees to adjust and align the display coordinate system of the image. Based on the preset resolution parameters, export the two-dimensional projected images from each viewpoint after calibrating the display coordinate system.
10. A photoacoustic vessel segmentation device based on the intersection of three-dimensional maximum intensity projections, characterized in that, The device includes a generation module, a construction module, and an output module; The generation module is used to acquire the three-dimensional photoacoustic image to be processed, and to project the three-dimensional photoacoustic image with maximum intensity along three mutually orthogonal axes to generate three corresponding two-dimensional projection images. The generation module is also used to determine the target blood vessel region on the two-dimensional projection image and generate three corresponding two-dimensional binary masks. The construction module is used to extend the two-dimensional binary mask to the three-dimensional voxel space along the corresponding orthogonal axes to construct three three-dimensional extended masks. The construction module is also used to perform three-dimensional spatial logical intersection processing on the three-dimensional extended mask to construct a three-dimensional blood vessel segmentation mask; The output module is used to extract the vascular structure from the three-dimensional photoacoustic image based on the three-dimensional vascular segmentation mask and output a preliminary segmentation image; The output module is also used to perform two-stage background purification on the preliminary segmented image and output a purified three-dimensional blood vessel segmentation image; wherein, the two-stage background purification includes isolated noise removal processing based on spatial local connectivity features and background truncation processing based on global signal intensity.