A target point driven arbitrary cross-sectional CT image visualization method and system

Abstract

The application discloses a target point driven arbitrary cross-section CT image visualization method and system, and the method comprises the following steps: receiving a selected point coordinate, using the speed and acceleration of a historical trajectory to predict a depth mapping to obtain an initial target point; taking the initial target point as a starting mass point, calculating a virtual gravity resultant force based on a gray histogram multi-peak detection and moving the mass point, determining an anchoring target point through oscillation suppression; taking the anchoring target point as a cross-section origin, generating a cross-section according to a gradient direction histogram adaptive fusion user normal vector; constructing a virtual spring system based on a projection point offset to correct a rotation angular velocity, updating the cross-section posture; adopting a neighborhood prediction and segmented variable step length strategy to advance a sampling ray along the normal vector, and stopping when the accumulated gray difference reaches a threshold value; finally, generating a cross-section image through median filtering, adaptive thickness Gaussian weighting and contrast stretching resampling. Precise positioning, adaptive cross-section orientation, efficient sampling and high-quality imaging are realized, and the interactive efficiency is significantly improved.

Description

Technical Field

[0001] This invention belongs to the field of medical image processing technology, and in particular relates to a target point-driven arbitrary cross-sectional CT image visualization method and system. Background Technology

[0002] In medical imaging diagnosis, computed tomography (CT) images provide three-dimensional information about the internal structures of the human body. When reviewing three-dimensional volume data, doctors often need to observe cross-sectional images of arbitrary oblique sections to examine the relationship between lesions and surrounding tissues from different angles. Traditional methods for generating arbitrary cross-sections usually require users to manually adjust the position and orientation of the cross-section, making the interactive process cumbersome and difficult to accurately locate the anatomical structures of interest.

[0003] Existing technologies include several methods for cross-section generation based on target point localization. Users select a point of interest by clicking on a location in a 2D display window, and the system maps this point to 3D volume data, generating a cross-section using this point as the origin. However, due to the lack of depth information in 2D clicking, the initial target point obtained is often inaccurate and deviates from the anatomical structure the user is truly interested in. Furthermore, the setting of the cross-section orientation lacks adaptability to local anatomical structures, requiring users to repeatedly adjust the normal vector to obtain a satisfactory cross-section. During sampling, traditional methods use equal-step progression for all pixels, resulting in low efficiency. Therefore, there is an urgent need for a visualization method that can accurately locate target points, adaptively determine the cross-section orientation, efficiently sample, and generate high-quality cross-section images. Summary of the Invention

[0004] This invention provides a target-point driven arbitrary section CT image visualization method and system to solve the technical problems of inaccurate target point positioning, poor adaptability to section direction and low sampling efficiency in the prior art.

[0005] In a first aspect, the present invention provides a target-point driven method for visualizing arbitrary cross-sectional CT images, comprising: The target medical image sequence is acquired, and three-dimensional volume data is constructed based on the spatial location information and grayscale value of each frame of the target medical image. The three-dimensional volume data is then uploaded to the video memory of the graphics processor to form a volume data resource that can be shared by any cross-sectional view. Receive the selected point coordinates input by the user, and map the selected point coordinates to the three-dimensional volume data according to the preset coordinate mapping strategy to obtain the initial target point; Using the initial target point as the starting mass, the gray value of each voxel in the three-dimensional volume data is used as the virtual mass. The virtual gravitational resultant force of all voxels in the preset spherical neighborhood centered on the current position is calculated. The starting mass is moved along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold. The final stopping position is defined as the anchoring target point. Using the anchored target point as the origin of the tangent, generate the current arbitrary cross section based on the initial normal vector and the in-plane direction vector; The rotational angular velocity input by the user is obtained. An inertia coefficient is generated based on the offset vector between the projection point of the anchored target point on the current arbitrary cross-section and the geometric center of the current arbitrary cross-section. The rotational angular velocity is corrected using the inertia coefficient, and the spatial attitude of the current arbitrary cross-section is updated to obtain the current initial arbitrary end face. Define a sampling ray for each pixel in the current initial arbitrary end face, and advance point by point along each sampling ray from the starting point. Calculate the absolute difference in gray level between the current sampling point and the previous sampling point and accumulate the difference. Stop sampling the sampling ray when the accumulated value reaches a preset threshold. The volume data resource is resampled based on the sampling results of all sampled rays to generate an arbitrary cross-section of the current target, and the cross-sectional image corresponding to the arbitrary cross-section of the current target is displayed in the arbitrary cross-sectional view.

[0006] In a second aspect, the present invention provides a target-point driven arbitrary cross-sectional CT image visualization system, comprising: The acquisition module is configured to acquire the target medical image sequence, construct three-dimensional volume data based on the spatial location information and grayscale value of each frame of the target medical image, and upload the three-dimensional volume data to the graphics processor's video memory to form a volume data resource that can be shared by any cross-sectional view. The mapping module is configured to receive the selected point coordinates input by the user and map the selected point coordinates to the three-dimensional volume data according to a preset coordinate mapping strategy to obtain the initial target point; The calculation module is configured to take the initial target point as the starting mass, use the gray value of each voxel in the three-dimensional volume data as the virtual mass, calculate the virtual gravitational resultant force of all voxels in the preset spherical neighborhood centered on the current position on the starting mass, and move the starting mass along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold, and define the final stopping position as the anchoring target point. The generation module is configured to use the anchored target point as the origin of the tangent and generate any current cross section based on the initial normal vector and the in-plane direction vector. The update module is configured to obtain the rotational angular velocity input by the user, generate an inertia coefficient based on the offset vector between the projection point of the anchored target point on the current arbitrary cross-section and the geometric center of the current arbitrary cross-section, use the inertia coefficient to correct the rotational angular velocity, and update the spatial attitude of the current arbitrary cross-section to obtain the current initial arbitrary end face. The sampling module is configured to define a sampling ray for each pixel in the current initial arbitrary end face, advance along each sampling ray from the starting point point, calculate the absolute difference in gray level between the current sampling point and the previous sampling point in turn and accumulate it, and stop sampling the sampling ray when the accumulated value reaches a preset threshold. The output module is configured to resample the volume data resource based on the sampling results of all sampled rays, generate an arbitrary cross-section of the current target, and display the cross-sectional image corresponding to the arbitrary cross-section of the current target in the arbitrary cross-sectional view.

[0007] Thirdly, an electronic device is provided, comprising: at least one processor, and a memory communicatively connected to the at least one processor, wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to perform the steps of the target point-driven arbitrary section CT image visualization method of any embodiment of the present invention.

[0008] Fourthly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the program instructions are executed by a processor, the processor performs the steps of the target point-driven arbitrary section CT image visualization method of any embodiment of the present invention.

[0009] This application presents a target-point-driven arbitrary cross-sectional CT image visualization method and system. First, it constructs 3D volume data and uploads it to GPU memory. Then, it receives the coordinates of selected points and uses the velocity and acceleration of historical trajectories to predict depth mapping to obtain an initial target point. Using the initial target point as the starting mass, it calculates the virtual gravitational resultant force based on multi-peak detection of the grayscale histogram and moves the mass, determining the anchoring target point through oscillation suppression. Using the anchoring target point as the origin of the cross-section, it adaptively fuses the user's normal vector based on the gradient direction histogram to generate the cross-section. Based on the projection point offset, it constructs a virtual spring system to correct the rotational angular velocity and update the cross-section attitude. It employs a neighborhood prediction and piecewise variable step-size strategy to advance the sampling ray along the normal vector, stopping when the accumulated grayscale difference reaches a threshold. Finally, it generates the cross-sectional image through median filtering, adaptive thickness Gaussian weighting, and contrast stretching resampling. This achieves precise positioning, adaptive cross-sectional orientation, efficient sampling, and high-quality imaging, significantly improving interaction efficiency. Attached Figure Description

[0010] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0011] Figure 1 A flowchart illustrating a target-point-driven arbitrary section CT image visualization method according to an embodiment of the present invention; Figure 2 This is a structural block diagram of a target-point driven arbitrary section CT image visualization system provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. Detailed Implementation

[0012] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0013] Please see Figure 1 The flowchart illustrates a target-point driven arbitrary cross-sectional CT image visualization method according to this application.

[0014] like Figure 1 As shown, the target-point driven arbitrary section CT image visualization method specifically includes the following steps: Step S101: Obtain the target medical image sequence, construct three-dimensional volume data based on the spatial location information and grayscale value of each frame of the target medical image, and upload the three-dimensional volume data to the graphics processor's video memory to form a volume data resource that can be shared by any cross-sectional view.

[0015] In this step, the DICOM format CT medical image sequence is first read and sorted in ascending order by slice position (e.g., InstanceNumber). For each frame, the pixel spacing, slice thickness, image size, and grayscale scaling parameters are extracted, and the Henle unit (HU) of each pixel is calculated. A three-dimensional array of width × height × number of frames is allocated in CPU memory, and the converted HU values ​​are stored in the corresponding voxel positions frame by frame to construct the three-dimensional volume data.

[0016] Then, a 3D texture object is created in GPU memory using a graphics API (such as OpenGL or CUDA), with the texture format set to 16-bit or 32-bit floating-point to maintain precision. The volume data from the CPU is transferred to this 3D texture all at once, with the texture filtering method set to trilinear interpolation and the addressing mode set to clamp to boundary. This texture resource is shared by all arbitrary section views; each view only needs to bind to the same texture to sample volume data in real time, eliminating the need for repeated storage. When loading new data, the existing texture is released before reconstruction.

[0017] Step S102: Receive the selected point coordinates input by the user, and map the selected point coordinates to the three-dimensional volume data according to the preset coordinate mapping strategy to obtain the initial target point.

[0018] In this step, the input time of the selected point coordinates is obtained, and the coordinates of multiple consecutive historical selected points before the input time are recorded to form a coordinate sequence; The user's line of sight is fitted on the display window plane according to the coordinate sequence, and the instantaneous velocity direction and instantaneous acceleration direction of the movement trajectory at the current moment are calculated. The instantaneous velocity direction and instantaneous acceleration direction are both two-dimensional direction vectors in the display window plane. Starting from the selected point coordinates, a ray is projected into the display window plane along the instantaneous velocity direction, and the ray is extended from the display window into the three-dimensional scene to form the first spatial ray. The first intersection point of the first spatial ray and the spatial bounding box of the three-dimensional volume data is calculated as the first intersection point, and the distance between the first intersection point and the observation point is used as the first depth candidate value. The observation point is the location of the virtual camera in the three-dimensional scene. Starting from the selected point coordinates, another ray is projected into the display window plane along the instantaneous acceleration direction, and the other ray is extended from the display window into the three-dimensional scene to form a second spatial ray. The first intersection point of the second spatial ray and the spatial bounding box is calculated as the second intersection point, and the distance between the second intersection point and the observation point is used as the second depth candidate value. Calculate the angle between the instantaneous velocity direction and the instantaneous acceleration direction, and determine whether the angle is less than a preset angle threshold; If the included angle is less than the preset angle threshold, the weighted average of the first depth candidate value and the second depth candidate value will be used as the initial depth value of the initial target point. If the included angle is not less than the preset angle threshold, then the first depth candidate value is used as the initial depth value of the initial target point; The initial target point is constructed by using the selected point coordinates as two-dimensional plane coordinates and the initial depth value as spatial coordinates.

[0019] In one specific embodiment, the system receives the coordinates of a selected point entered by the user via a mouse or touchscreen on the display window of any cross-sectional view, and maps the selected point coordinates to the three-dimensional volume data according to a preset coordinate mapping strategy to obtain the initial target point.

[0020] First, the system records the moment the user selects a point and retrieves the coordinates of N consecutive historical points selected before that moment (e.g., N=5), forming a coordinate sequence. These historical point coordinates record the trajectory of the mouse movement on the display window before the user clicks the current point.

[0021] Then, the system fits the user's line of sight movement trajectory on the display window plane based on the coordinate sequence. Common fitting methods include least squares or Kalman filtering. By calculating the first and second derivatives of the trajectory, the instantaneous velocity direction and instantaneous acceleration direction of the movement trajectory at the current moment are calculated, respectively. The instantaneous velocity direction reflects the main trend of the user's current mouse movement, while the instantaneous acceleration direction reflects the changing trend of the user's movement intention (e.g., sudden turning or acceleration). Both are two-dimensional direction vectors within the display window plane.

[0022] Example: Suppose the coordinates of the user's five consecutive historical selected points are (100,100), (105,102), (112,105), (120,108), and (130,112). Through difference calculation, the instantaneous velocity direction is approximately (1,0.4) (after normalization), and the instantaneous acceleration direction is approximately (0.2,0.05) (after normalization), indicating that the user is moving rapidly to the upper right with a gradual change in speed.

[0023] Next, the system projects a ray along the instantaneous velocity direction onto the display window plane, starting from the current selected point coordinates, and extends this ray from the 2D display window into the 3D scene, forming a first spatial ray. The first intersection point of this first spatial ray with the spatial bounding box of the 3D volume data (i.e., the point where the ray first enters the volume data bounding box from the viewpoint) is calculated and recorded as the first intersection point. Simultaneously, the spatial straight-line distance between the viewpoint (i.e., the location of the virtual camera in the 3D scene) and the first intersection point is obtained, and this distance is used as the first depth candidate value. This depth candidate value reflects the approximate depth of the tissue structure pointed to by the user's current velocity direction.

[0024] Similarly, the system uses the current selected point coordinates as the starting point, projects another ray along the instantaneous acceleration direction into the display window plane, and extends it into the 3D scene to form a second spatial ray. It calculates the first intersection point (the second intersection point) with the spatial bounding box and the distance from that intersection point to the observation point, using this as the second depth candidate value. This depth candidate value reflects the depth indicated by the change in the user's movement intention.

[0025] In practical applications, to improve computational efficiency, ray-bounding box intersection algorithms (such as the slab method) can be used to quickly find the intersection points.

[0026] The system then calculates the angle between the instantaneous velocity direction and the instantaneous acceleration direction. This angle reflects the "smoothness" of the user's mouse movement: the smaller the angle, the more stable the user's movement direction and the clearer the intention; the larger the angle, the more likely the user is turning or accelerating, and the reliability of the velocity direction decreases.

[0027] The specific judgment logic is as follows: If the included angle is less than a preset angle threshold (e.g., 30°), it indicates that the velocity direction is consistent with the acceleration direction, and the user's depth intention is relatively certain. In this case, the weighted average of the first and second depth candidate values ​​is used as the initial depth value of the initial target point. The weighting coefficient can be set according to the actual situation, for example, the weight of the first depth candidate value is 0.7 and the weight of the second depth candidate value is 0.3, to highlight the dominant role of the velocity direction. In a preferred example, the weighting coefficient is inversely proportional to the size of the included angle: the smaller the included angle, the greater the weight of the first depth candidate value.

[0028] If the included angle is not less than the preset angle threshold, it indicates that the user's movement trajectory has changed drastically and the acceleration direction may be unreliable. In this case, the first depth candidate value is directly used as the initial depth value to avoid introducing errors.

[0029] Finally, the system uses the user-inputted coordinates of the selected point as the two-dimensional plane coordinates of the initial target point (i.e., the projected plane coordinates along the viewing direction in the three-dimensional volume data), and uses the initial depth value calculated above as the component along the viewing direction in the spatial coordinates. A three-dimensional spatial point is then constructed in the viewing coordinate system, thus obtaining the initial target point. If the volume data coordinate system is inconsistent with the viewing coordinate system, corresponding coordinate transformations are required.

[0030] In summary, by utilizing the historical trajectory information of continuous point selection operations by the user, depth is dynamically inferred from 2D point selection, effectively solving the technical problems of traditional single-point mapping's inability to determine depth and the arbitrary placement of target points on the volume data surface or fixed depth plane. Specifically: By statistically analyzing the user's movement trajectory over a short period of time and calculating instantaneous velocity and acceleration, the system can adaptively capture the user's depth intent. When the user slowly moves the mouse to point at a certain anatomical structure, the velocity direction is stable, and the accuracy of the depth candidate value output by the system is high. When the user clicks quickly or the trajectory changes abruptly, the system automatically switches to the first depth candidate value to avoid errors introduced by acceleration noise. By introducing angle judgment and weighted averaging mechanisms, the system can stably output reasonable depth values ​​under different operating habits (fast selection, slow aiming, and curve movement), overcoming the ambiguity problem of multiple candidate depths in the overlapping area of ​​the line of sight and volume data in the single ray method.

[0031] Step S103: Taking the initial target point as the starting mass, the gray value of each voxel in the three-dimensional volume data is used as the virtual mass. The virtual gravitational resultant force of all voxels in the preset spherical neighborhood centered on the current position is calculated. The starting mass is moved along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold. The final stopping position is defined as the anchoring target point.

[0032] In this step, the particle at the current iteration position is taken as the current particle, where the current particle is the starting particle in the first iteration; With the current position of the current particle as the center point, construct a spherical neighborhood, and the radius of the spherical neighborhood is the preset initial radius; Obtain the gray values ​​of all voxels within the spherical neighborhood, and statistically analyze the distribution of gray values ​​to determine the first gray value with the highest frequency and the second gray value with the second highest frequency. Calculate the absolute difference between the first gray value and the second gray value, and determine whether the absolute difference is greater than the preset separation threshold; If the absolute difference is greater than the preset separation threshold, it is determined that there are at least two different tissue types in the spherical neighborhood, and the gray centroid corresponding to each tissue type is calculated respectively. Among them, voxels with gray values ​​within the first preset gray value range are classified as the first tissue type, and voxels with gray values ​​within the second preset gray value range are classified as the second tissue type. Calculate the first gravitational vector pointing from the current particle to the gray-scale centroid of the first tissue type, and the second gravitational vector pointing from the current particle to the gray-scale centroid of the second tissue type. The magnitude of the first gravitational vector is proportional to the proportion of the number of voxels contained in the first tissue type to the total number of voxels in the spherical neighborhood, and the magnitude of the second gravitational vector is proportional to the proportion of the number of voxels contained in the second tissue type to the total number of voxels in the spherical neighborhood. By superimposing the first gravitational vector and the second gravitational vector, the virtual gravitational resultant force on the current particle is obtained; If the absolute difference is not greater than the preset separation threshold, the spherical neighborhood is determined to be a single tissue type. The gray value of each voxel in the spherical neighborhood is normalized to the range of zero to one and used as the virtual mass of the voxel. Then, the gravitational vector of each voxel to the center point is calculated. All gravitational vectors are summed to obtain the virtual gravitational resultant force. The direction of the gravitational vector of each voxel is the direction from the voxel to the center point. The magnitude of the gravitational vector is directly proportional to the virtual mass of the voxel and inversely proportional to the square of the distance from the voxel to the center point.

[0033] It should be noted that for the current particle at the current iteration position, the current particle in the first iteration is the starting particle, and the current position of the current particle and the current virtual gravitational resultant force at the current position are obtained; Determine whether the magnitude of the current virtual gravitational resultant force is less than the preset gravitational threshold. If the force is less than the preset gravity threshold, the movement will stop immediately, and the current position of the particle will be defined as the anchor target point. If it is not less than the preset gravity threshold, then obtain the historical position before the last movement and calculate the historical displacement vector from the historical position to the current particle. Calculate the angle between the historical displacement vector and the current virtual gravitational resultant force direction, and determine whether the angle is greater than the preset oscillation angle threshold in each of the preset iterations. If the included angle is greater than the preset oscillation angle threshold in consecutive preset iterations, it is determined that the current particle is trapped in a local oscillation state, and an oscillation suppression operation is performed to obtain a new particle position. Otherwise, the direction of the current virtual gravitational resultant force is used as the moving direction, and the current particle is moved with the current moving step size to obtain a new particle position. The current moving step size is proportional to the magnitude of the current virtual gravitational resultant force. The new particle position is used as the current particle position in the next iteration. The virtual gravitational resultant force is recalculated until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold. At this point, the position of the current particle is defined as the anchor target point.

[0034] Specifically, performing oscillation suppression operations yields new particle positions including: Record the direction of the oscillation axis, which is the average direction of multiple consecutive historical displacement vectors; The current virtual gravitational force is decomposed into a first component along the oscillation axis and a second component perpendicular to the oscillation axis. Multiply the first component by the preset damping coefficient to obtain the suppressed first component; The suppressed first component and the second component are recombined to obtain the suppressed virtual gravitational resultant force; Using the direction of the suppressed virtual gravitational net force as the direction of movement, and taking half of the current movement step as the new movement step, move the current particle along the direction of movement to obtain the new position of the particle.

[0035] In one specific embodiment, the initial target point obtained in step S102 is used as the starting particle. By simulating the movement of the particle under virtual gravity in the gray field of three-dimensional volume data, the particle is eventually stabilized in a local gray-level extreme region (usually the center of high-density tissue), and this position is defined as the anchoring target point. This process is divided into two sub-steps: virtual gravitational resultant force calculation and iterative movement of the particle.

[0036] I. Calculation of Virtual Gravitational Resultant Force The system uses the particle at the current iteration position as the current particle. In the first iteration, the current particle is the initial target point obtained in step S102.

[0037] A spherical neighborhood is constructed with the current position of the current particle as the center point. The radius of this spherical neighborhood is a preset initial radius, such as 3 voxel spacings. The selection of this radius requires a trade-off between computational efficiency and local feature capture capability: a radius that is too small is prone to getting trapped in local noise, while a radius that is too large will introduce too many distant voxels, causing a smooth transition in the gravitational direction. In a preferred example, the initial radius is adaptively set according to the spatial resolution of the volume data, for example, 2 to 3 times the average voxel spacing.

[0038] Obtain the grayscale values ​​of all voxels within the spherical neighborhood and statistically analyze their distribution. A specific statistical method can be using a grayscale histogram, dividing the grayscale range into several intervals (e.g., 256 intervals), and counting the number of voxels in each interval. Then, determine the most frequent grayscale value (i.e., the peak value of the histogram) and the second most frequent grayscale value (i.e., the second peak value of the histogram).

[0039] The absolute difference between the first and second grayscale values ​​is calculated, and it is determined whether this absolute difference is greater than a preset separation threshold. The separation threshold is related to the grayscale resolution of the CT image. In one example, if the grayscale value range is 0~4095 (12-bit CT image), the separation threshold can be set to 100~200. This determination is used to identify whether the current spherical neighborhood contains multiple tissue types: if the difference is greater than the threshold, it indicates the existence of two clearly separated grayscale peaks; if the difference is not greater than the threshold, it indicates that the grayscale distribution is concentrated, possibly representing a single tissue or noise.

[0040] Scenario 1: Multiple organizational types exist If the absolute difference is greater than a preset separation threshold, it is determined that at least two different tissue types (e.g., bone and muscle, tumor and normal tissue) exist within the spherical neighborhood. The system calculates the gray-level centroid corresponding to each tissue type. Specifically: Voxels whose gray values ​​fall within a preset range of the first gray value (e.g., the first gray value ± threshold δ1) are classified as the first tissue type. These voxels typically correspond to tissues with higher density (such as bones).

[0041] Voxels whose gray values ​​fall within the preset range of the second gray value (e.g., the second gray value ± threshold δ2) are classified as the second tissue type. These voxels typically correspond to tissues with lower density (e.g., soft tissue).

[0042] For each tissue type, the gray centroid of that tissue type is calculated by weighting the spatial coordinates of all voxels contained therein according to their gray values. The weighting method can use gray values ​​as weights to make high-gray-level voxels contribute more to the centroid.

[0043] Then, the first gravitational vector pointing from the current particle to the gray-scale centroid of the first tissue type, and the second gravitational vector pointing from the current particle to the gray-scale centroid of the second tissue type are calculated respectively. The magnitude (modulus) of the gravitational vector is proportional to the proportion of the number of voxels contained in the corresponding tissue type to the total number of voxels in the spherical neighborhood. For example, if the first tissue type contains 100 voxels and the total number of voxels in the spherical neighborhood is 500, the proportionality coefficient is 0.2. The larger this proportion, the more dominant the tissue is locally, and the stronger its gravitational pull on the particle.

[0044] The first and second gravitational vectors are superimposed to obtain the virtual resultant gravitational force acting on the current particle. The direction of this resultant force points towards the "comprehensive centroid" of the local multi-organism grayscale distribution, which helps the particle move from the organization boundary towards the interior of the main organization.

[0045] Scenario 2: Single Organization Type If the absolute difference is not greater than the preset separation threshold, the spherical neighborhood is determined to be a single tissue type (or a region with a gradual change in grayscale). The system normalizes the grayscale value of each voxel in the spherical neighborhood to the range of zero to one, and uses this as the virtual quality of that voxel. The normalization formula is: Virtual quality = (Current grayscale value - Minimum grayscale value in the neighborhood) / (Maximum grayscale value in the neighborhood - Minimum grayscale value in the neighborhood), so that high grayscale voxels have greater virtual quality.

[0046] For each voxel within the spherical neighborhood, calculate the gravitational vector of that voxel to the center point (i.e., the current position of the particle). The direction of the gravitational vector is from the voxel towards the center point. The magnitude of the gravitational vector is directly proportional to the virtual mass of the voxel and inversely proportional to the square of the distance from the voxel to the center point. This design simulates the law of universal gravitation: the closer the voxel is to the center point and the greater its mass, the stronger its attraction to the particle.

[0047] The virtual gravitational force acting on the current particle is obtained by summing the gravitational vectors of all voxels. This resultant force points towards the centroid of the grayscale distribution in the neighborhood, causing the particle to naturally move towards the center of the high grayscale region.

[0048] II. Iterative Movement and Oscillation Suppression of Particles After obtaining the virtual gravitational net force of the current particle, the system determines whether the particle needs to be moved and how to move it.

[0049] First, determine whether the magnitude of the current virtual gravitational net force is less than the preset gravitational threshold. This threshold is a very small positive number (e.g., 0.01). When the net force is close to zero, it means that the particle has reached a local gravitational equilibrium position (i.e., a local extreme point in the grayscale distribution). At this point, stop moving directly and define the current position of the particle as the anchor target point.

[0050] If the resultant force modulus is not less than the threshold, the particle needs to continue moving. The system obtains the historical position before the last movement and calculates the historical displacement vector from the historical position to the current particle. This displacement vector reflects the direction and distance of the last movement.

[0051] Calculate the angle between the historical displacement vector and the current direction of the virtual gravitational resultant force. This angle reflects the consistency between the particle's direction of movement and the direction of gravity: if the angle is small, it indicates that the particle is moving steadily along the direction of the resultant force; if the angle is close to 90° or larger, it indicates that the particle may be oscillating back and forth.

[0052] The system determines whether the included angle is greater than a preset oscillation angle threshold (e.g., 60°) in each of the preset number of iterations (e.g., 3 consecutive iterations). If this condition is met, and the magnitude of the virtual gravitational force is observed to alternately increase and decrease in adjacent iterations (i.e., exhibit fluctuations), then the current particle is determined to be in a local oscillation state. This oscillation usually occurs in a flat region near the target point or around an extreme point caused by discrete sampling.

[0053] Oscillation suppression operation: Record the oscillation axis direction: The oscillation axis direction is the average direction (normalized) of multiple consecutive historical displacement vectors. This direction represents the axis of repeated oscillation of the particle.

[0054] The current virtual gravitational force is decomposed into a first component along the oscillation axis and a second component perpendicular to the oscillation axis. This can be achieved through vector dot product projection.

[0055] Multiplying the first component by a preset damping coefficient (e.g., 0.3~0.5) yields the suppressed first component. A damping coefficient less than 1 is used to weaken the gravitational force along the oscillation axis, thereby reducing the oscillation amplitude.

[0056] The suppressed first component and the second component are recombined to obtain the suppressed virtual gravitational resultant force.

[0057] Using the direction of the suppressed virtual gravitational net force as the direction of movement, and taking half of the current movement step as the new movement step, move the current particle along this direction to obtain the new particle position.

[0058] Normal movement when no oscillation is detected: If the included angle does not meet the oscillation condition (i.e., at least once the included angle is not greater than the threshold) in consecutive preset iterations, then the direction of the current virtual gravitational resultant force is directly used as the movement direction, and the current particle is moved with the current movement step size to obtain the new particle position. The current movement step size is proportional to the magnitude of the current virtual gravitational resultant force. A larger magnitude indicates that the particle is far from the equilibrium point, and a larger movement step should be taken for rapid convergence; a smaller magnitude indicates that the particle is close to the equilibrium point, and a smaller movement step should be taken to improve accuracy. For example, the step size = coefficient × resultant force magnitude, and the maximum step size is limited to no more than half the radius of the spherical neighborhood.

[0059] After obtaining the new particle position, use it as the current particle position for the next iteration, and repeat the above gravitational calculation and movement process. Continue iterating until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold. At this point, the position of the current particle is defined as the anchoring target point.

[0060] In a specific clinical application scenario, the user wants to locate the center of a liver tumor lesion. The initial target point may fall on the tumor margin or nearby normal tissue. Through the aforementioned gravitational movement process, the particle is attracted by the high-grayscale region (the tumor is usually slightly higher or lower in density than normal liver tissue, depending on the enhancement method) and gradually moves towards the tumor center. Even if there is a partial volume effect or noise at the tumor boundary causing slight oscillations in the particle, the oscillation suppression mechanism can quickly smooth out the fluctuations, allowing the particle to converge stably to the core of the tumor. The final anchored target point is the origin of the section desired by the physician.

[0061] In summary, multi-peak detection using grayscale histograms distinguishes different tissue types, and the weighted gravity of the centroid of each tissue is calculated separately, enabling particles to accurately move into the target tissue rather than remaining at the tissue boundary. Compared to traditional gradient descent or simple centroid methods, the localization accuracy is significantly improved. For single tissue regions with uniform grayscale distribution (such as liver parenchyma), an inverse square distance gravity model is used to move particles towards the grayscale centroid. For multi-tissue boundary regions, multi-target gravity synthesis is used, avoiding the directional ambiguity of a single gravity at the boundary. Regardless of whether the target anatomical structure is homogeneous or heterogeneous, this method can reliably converge.

[0062] Step S104: Using the anchored target point as the origin of the tangent, generate the current arbitrary cross-section based on the initial normal vector and the in-plane direction vector.

[0063] In this step, any cross section contains a two-dimensional pixel grid determined by pixel width, pixel height and sampling interval. Each pixel in the two-dimensional pixel grid corresponds to a three-dimensional spatial point in the three-dimensional volume data.

[0064] Specifically, the local neighborhood at the anchor target point is obtained. The local neighborhood is a cubic region with a preset side length centered on the anchor target point. Calculate the gray-level gradient of each voxel in the local neighborhood to obtain the gradient direction vector of each voxel, and count the direction distribution of all gradient direction vectors to construct a gradient direction histogram. The horizontal axis of the gradient direction histogram is the direction angle interval, and the vertical axis is the number of gradient direction vectors falling within the direction angle interval. The principal gradient direction and the secondary gradient direction are identified from the gradient direction histogram. The principal gradient direction is the direction corresponding to the direction angle interval that appears most frequently, and the secondary gradient direction is the direction that appears most frequently in the direction angle interval that is orthogonal to the principal gradient direction. The subgradient direction is used as the candidate direction for the cross-sectional normal vector, and the principal gradient direction is used as the candidate direction for the in-plane direction vector. Calculate the target angle between the candidate direction and the initial normal vector input by the user; If the target angle is less than the preset first angle threshold, the candidate direction will be used as the normal vector of any current cross section. If the target angle is greater than the preset second angle threshold, the initial normal vector is used as the normal vector of any current cross section, where the second angle threshold is greater than the first angle threshold. If the included angle is between the first angle threshold and the second angle threshold, then the weighted average direction of the candidate direction and the initial normal vector is calculated as the normal vector of the current arbitrary section. Determine the target direction that is simultaneously perpendicular to the normal vector of the current arbitrary cross section and the preset global upward direction vector. Take the target direction as the in-plane direction vector and normalize the in-plane direction vector to obtain the normalized in-plane direction vector. The global upward direction vector is a unit vector pointing vertically upward in three-dimensional space. Using the anchored target point as the origin of the tangent, and based on the normal vector of the current arbitrary cross-section and the normalized in-plane direction vector, combined with the size of the two-dimensional pixel grid and the sampling interval, the coordinates of the three-dimensional spatial point corresponding to each pixel are calculated to generate the current arbitrary cross-section.

[0065] In one specific embodiment, firstly, a cubic region is defined as a local neighborhood centered on the anchor target point. The side length of the cube is preset to a fixed value based on the spatial resolution of the volume data, such as 10 mm or 20 voxel spacing, to ensure that the main structural features of the local tissue can be captured while avoiding the introduction of irrelevant information from afar.

[0066] For each voxel within the cubic region, the system calculates its gray-level gradient. The gray-level gradient reflects the degree and direction of change of gray values ​​around the voxel in space. Specifically, for each voxel, the rate of change of gray value (i.e., partial derivatives) along the X, Y, and Z axes is calculated, using either the central difference method or the Sobel convolution operator. The rates of change in the three directions are combined into a three-dimensional vector, which is the gradient direction vector of that voxel. The gradient direction points in the direction of the fastest increase in gray value. In CT images, the boundaries between different tissue types (such as bone and soft tissue, soft tissue and air) typically have larger gradient amplitudes, and the gradient direction is perpendicular to the tangent plane of the interface.

[0067] Since gradient directions are directions in three-dimensional space, they are difficult to count directly by angle intervals. Therefore, the system uses spherical coordinates to represent each gradient direction, uniquely identifying a direction using azimuth (0°~360°) and pitch (-90°~90°). To construct the histogram, the azimuth and pitch angles are discretized at certain intervals, for example, 15° intervals for both azimuth and pitch, forming a series of solid angle grids. For each voxel in a local neighborhood, its gradient direction is assigned to the corresponding solid angle grid only if its gradient magnitude exceeds a preset threshold (to exclude noise directions in flat areas). The number of voxels in each grid (i.e., frequency) is then counted. This yields the gradient direction histogram, with the horizontal axis representing the solid angle intervals and the vertical axis representing the frequency.

[0068] The system identifies the solid angle interval with the highest frequency from the histogram, and the direction represented by this interval is called the principal gradient direction. The principal gradient direction reflects the direction in which gray-level changes are most dominant in a local region, and is usually perpendicular to the main extension direction of the tissue structure (for example, for a blood vessel, the principal gradient direction is perpendicular to the vessel wall and points towards the lumen or the outside of the vessel). Next, the system finds the direction with the highest frequency among the solid angle intervals orthogonal to the principal gradient direction (i.e., the included angle is approximately 90°), and calls it the secondary gradient direction. The specific method is as follows: on the great circle of a sphere corresponding to the principal gradient direction (the circle formed by all directions perpendicular to the principal gradient direction), the frequency of each solid angle interval falling near this great circle is counted, and the direction with the highest frequency is taken as the secondary gradient direction. The secondary gradient direction is often parallel to the main extension direction of the tissue structure (e.g., along the long axis of a blood vessel or the longitudinal axis of a bone).

[0069] The system uses the subgradient direction as a candidate direction for the cross-sectional normal vector because this direction is parallel to the main extension direction of the tissue structure. A cross-section generated using this direction as the normal vector will be perpendicular to the main extension direction of the tissue (i.e., cut along the tissue's cross-section). Simultaneously, the main gradient direction is used as a candidate direction for the in-plane direction vector. This direction is perpendicular to the main extension direction of the tissue structure and is suitable as the in-plane direction, enabling the cross-sectional image to display the structure of the tissue's cross-section.

[0070] Users input an initial normal vector through the interactive interface, expressing the doctor's subjective expectation of the cross-sectional orientation (e.g., hoping to observe a standard coronal plane, sagittal plane, or a specific oblique section). The system calculates the target angle between the candidate direction (subgradient direction) and the initial normal vector, which is the minimum angle between the two direction vectors in three-dimensional space, ranging from 0° to 180°. This angle reflects the degree of difference between the natural orientation of the local anatomical structure and the user's desired direction.

[0071] The system presets two angle thresholds: a first angle threshold (a smaller value, e.g., 30°) and a second angle threshold (a larger value, e.g., 60°), with the first threshold being smaller than the second threshold. Based on the size of the target angle, three different decision strategies are used to determine the final cross-sectional normal vector: If the target angle is less than the first angle threshold, it means that the user's initial normal vector is very close to the direction of the local tissue structure. In this case, the anatomical structure information is used first, and the candidate direction (subgradient direction) is directly used as the final normal vector of any current cross-section. This allows the cross-sectional image to fit the tissue's own cross-section to the greatest extent, which is beneficial for showing the hierarchical relationship between structures.

[0072] If the target angle is greater than the second angle threshold, it means that the direction the user wants to observe is significantly different from the local anatomical structure (for example, to cut obliquely to expose a lesion that does not follow the normal direction). In this case, respect the user's interaction intention and use the initial normal vector input by the user as the final normal vector.

[0073] If the target angle falls between the first and second angle thresholds, the system uses a weighted average to fuse the candidate directions and the initial normal vector. Specifically, the weight of the candidate direction is set to a value obtained by linearly mapping the target angle: weight = (second angle threshold - target angle) / (second angle threshold - first angle threshold), and the weight of the initial normal vector is 1 minus this weight. The two direction vectors are multiplied by their corresponding weights, and then the two weighted vectors are added together. Finally, the sum is normalized (i.e., scaled to a length of 1) to obtain the weighted average direction as the final normal vector. This method achieves a smooth transition between anatomical constraints and user intent, avoiding the abrupt changes that binary selection might bring.

[0074] After the final normal vector is determined, it is also necessary to determine the row direction (i.e., the in-plane direction vector) within the cross-sectional image to uniquely fix the rotational orientation of the cross-section. The system predefines a global upward direction vector, which is a unit vector pointing vertically upward in three-dimensional space (e.g., along the positive Y-axis in the world coordinate system, typically pointing towards the patient's head or above the scanning table). This vector provides a unified "upward" reference for all cross-sections, ensuring that the display orientation of the cross-section remains consistent across different interactive scenarios.

[0075] The system calculates a direction that is perpendicular to both the final normal vector and the global upward direction vector, and uses this direction as the in-plane direction vector. Geometrically, this direction is perpendicular to the normal vector (ensuring it lies within the cross-sectional plane) and also perpendicular to the global upward direction (making the row direction of the cross-sectional image as consistent as possible with the "horizontal" direction). If the final normal vector and the global upward direction vector are exactly parallel (e.g., the cross-section is horizontal), the above perpendicularity condition cannot uniquely determine the direction. In this case, the system can either use the initial in-plane direction vector input by the user or choose any horizontal direction perpendicular to the normal vector (e.g., the X-axis direction) as a substitute. After obtaining the in-plane direction vector, it is normalized (i.e., its length is adjusted to 1) to obtain the normalized in-plane direction vector.

[0076] The two-dimensional pixel grid of any cross-section has defined dimensions: pixel width W (unit: pixels), pixel height H (unit: pixels), and sampling interval d (unit: millimeters / pixel). The sampling interval is usually set to a value close to the original voxel spacing of the volume data to ensure that the resolution of the cross-sectional image matches that of the original data.

[0077] The system uses the anchor target point as the origin of the cross section, which is usually placed at the geometric center of the pixel grid, i.e., at column (W-1) / 2 and row (H-1) / 2. This way, the doctor's focus can be kept in the center of the field of vision during the interaction.

[0078] For each pixel in the grid, its position is uniquely determined by its row index i (0 to H-1) and column index j (0 to W-1). The 3D coordinates of the corresponding point in space are calculated according to the following steps: Calculate the offset distance of this pixel relative to the center of the grid geometry: Row direction offset = (i - (H - 1) / 2) × d; Column offset = (j - (W - 1) / 2) × d; Calculate the column direction vector (i.e., the direction perpendicular to the row direction within the cross-section): The column direction vector is chosen to be perpendicular to both the final normal vector and the normalized in-plane direction vector (i.e., the two are orthogonal). This can be directly obtained through geometric relationships: the column direction vector, the normal vector, and the in-plane direction vector form a right-handed orthogonal system.

[0079] Calculate the 3D spatial coordinates of a pixel: Pixel coordinates = Anchored target point coordinates + Row direction offset × Normalized in-plane direction vector + Column direction offset × Column direction vector; The above calculations are performed sequentially on all pixels, and the resulting 3D spatial coordinates are organized in row and column order to generate the current arbitrary cross-section. This cross-section has a definite position and orientation in 3D space, laying the foundation for subsequent sampling ray advancement and cross-section image generation.

[0080] Application Example: In a liver CT visualization, the doctor wants to observe the cross-section along the portal vein. The doctor clicks on a point on the main trunk of the portal vein, and the system anchors this point to the center of the vessel in step S103. In step S104, the system calculates a gradient direction histogram in the local neighborhood around the anchor point: the gray-level gradient around the portal vein mainly reflects the boundary between the vessel wall and the liver parenchyma; the principal gradient direction is perpendicular to the vessel wall, and the secondary gradient direction is parallel to the long axis of the vessel. The initial normal vector input by the doctor is approximately along the long axis of the portal vein (e.g., at a 20° angle to the horizontal plane). The calculated angle between the candidate direction (secondary gradient direction) and the initial normal vector is 12°, which is less than the first threshold of 30°. The system adopts the candidate direction as the final normal vector, i.e., the direction of the long axis of the vessel. Then, the system determines a direction that is simultaneously perpendicular to this normal vector and the global upward direction as the in-plane direction vector, keeping the cross-sectional direction horizontal. Finally, a 512×512 pixel cross-sectional grid with a sampling interval of 1 mm is generated centered on the anchor point to obtain an oblique section along the portal vein, clearly showing the portal vein and its branches. Doctors can obtain an ideal diagnostic cross-section without manually adjusting the angle.

[0081] Step S105: Obtain the rotational angular velocity input by the user, generate an inertia coefficient based on the offset vector between the projection point of the anchored target point on the current arbitrary cross-section and the geometric center of the current arbitrary cross-section, use the inertia coefficient to correct the rotational angular velocity, and update the spatial attitude of the current arbitrary cross-section to obtain the current initial arbitrary end face.

[0082] In this step, the vertical projection point of the anchoring target point on the plane of the current arbitrary section is calculated; Obtain the geometric center point of any current cross-section display window, calculate the offset vector from the geometric center point to the vertical projection point, and calculate the magnitude of the offset vector; Obtain the diagonal length of the display window and calculate the normalized offset distance; A one-dimensional virtual spring-mass system is constructed with the geometric center point as the origin and the direction pointed to by the offset vector as the positive direction. The vertical projection point is regarded as a mass point connected to the spring, the geometric center point is the fixed point at the other end of the spring, and the spring constant is proportional to the normalized offset distance. According to Hooke's Law, calculate the restoring force on the particle and, based on the preset virtual particle mass, calculate the restoring acceleration. Get the scalar value of the rotational angular velocity input by the user, and get the rotational axis direction input by the user. The rotational axis direction is perpendicular to the current arbitrary cross-section and passes through the geometric center of the display window. The recovery acceleration is converted into recovery angular velocity. The conversion method is as follows: the magnitude of the recovery angular velocity is equal to the recovery acceleration divided by the diagonal length of the display window, and the direction of the recovery angular velocity is consistent with the direction of the rotation axis. The corrected rotational angular velocity is obtained by weighted fusion of the recovery angular velocity and the user-input rotational angular velocity. Based on the corrected rotational angular velocity and rotation axis direction, calculate the rotation update amount of the normal vector and in-plane direction vector of the current arbitrary cross-section to obtain the new spatial attitude, and define the updated cross-section as the current initial arbitrary end face.

[0083] In one specific embodiment, the system first calculates the vertical projection point of the anchored target point onto the plane containing the current arbitrary cross-section. Since the current arbitrary cross-section is uniquely determined by the normal vector and the origin of the tangent plane, the foot of the perpendicular from the anchored target point (a point in three-dimensional space) to this plane is the vertical projection point, denoted as point P. This projection point represents the position of the target structure of interest to the user on the cross-sectional image.

[0084] Next, the geometric center point of any current cross-sectional display window is obtained, denoted as point C. The display window is usually a rectangular area, and its geometric center is the center pixel position of the window. The system calculates the offset vector from the geometric center point C to the vertical projection point P, denoted as vector V=PC. At the same time, the magnitude of this offset vector is calculated, denoted as |V|. This magnitude reflects the distance of the anchor point that the user is focusing on from the center of the field of view: the larger |V| is, the more serious the deviation; the smaller |V| is, the closer to the center.

[0085] To facilitate consistent processing across display windows of different sizes, the system obtains the diagonal length of the display window, denoted as L. The normalized offset distance ε = |V| / L is calculated. The value of ε ranges from 0 to 1, where 0 indicates that the projection of the anchor point falls exactly at the center of the window, and 1 indicates that the projection point is located at a corner of the window.

[0086] The system constructs a one-dimensional virtual spring-mass system with the geometric center point C as the origin and the direction pointed to by the offset vector V as the positive direction. In this system: The vertical projection point P is considered as a point mass connected to the end of the spring; The geometric center point C is considered as the fixed point (immobile point) at the other end of the spring. The spring constant k is proportional to the normalized offset distance ε, i.e., k = k0 × ε, where k0 is a preset proportionality constant. This design makes the spring pull stronger when the deviation is greater, creating an interactive feeling similar to "gravity".

[0087] According to Hooke's Law (the restoring force of a spring is proportional to its deformation), the magnitude of the restoring force F on point P is: F = k × |V| = k0 × ε × |V|. Since ε = |V| / L, therefore F = k0 × (|V|²) / L. The direction of the restoring force always points from point P to point C (i.e., opposite to the direction of the offset vector V). The function of the restoring force is to pull the point mass back to its equilibrium position (the center of the window).

[0088] To convert the restoring force into a kinematic quantity, the system presets a virtual mass m (for example, a value of 1, representing unit mass, for ease of calculation). According to Newton's second law, the magnitude of the restoring acceleration a of mass P is: a = F / m = (k0 × ε × |V|) / m. The direction of acceleration is also from P to C.

[0089] Users input rotation via mouse dragging, touch swiping, or peripherals (such as a trackball). The system captures the scalar value ω of the rotational angular velocity input by the user. user and the direction of rotation axis R axis Rotation axis direction R axis The geometric meaning of is: the axis perpendicular to any current cross-section and passing through the geometric center point of the display window. This is because in cross-section browsing, rotation usually occurs around the center of the view, making the cross-section appear to rotate around an axis passing through the center. The direction of the rotation axis is determined by the cross-section normal vector, i.e., R. axis Parallel to the cross-section normal vector (same or opposite, depending on the convention of rotation direction).

[0090] The restoring force and restoring acceleration in the system are linear kinematic quantities, while the rotational angular velocity is an angular kinematic quantity. Therefore, dimension matching and conversion are required. The system converts the restoring acceleration into the restoring angular velocity ω. recovery The conversion method is as follows: the magnitude of the recovery angular velocity is equal to the recovery acceleration a divided by the diagonal length L of the display window, i.e., |ω recovery |=a / L; The direction of the recovery angular velocity is relative to the rotation axis direction R axis Consistent. The physical meaning of this transformation is that it "maps" the linear recovery motion to a rotational motion around the center of the field of view, making the intuitiveness of the cross-sectional rotation proportional to the offset distance.

[0091] Weighted fusion: The system will recover the angular velocity ω recovery The rotational angular velocity ω input by the user user By performing a weighted summation, the corrected rotational angular velocity ω is obtained. final Weighted fusion uses a simple addition method: ,in This is a preset inertial fusion coefficient (typically 0.2~0.8, adjustable based on user interaction feel). After fusion, when the user actively rotates the cross-section, if the anchor point deviates from the center, the system will add a corrective bias, causing the cross-section to rotate slightly towards the center; when the user releases the mouse (ω... user Even when the center angle is 0, the system will continue to rotate based on the recovery angular velocity until the anchor point returns to the vicinity of the center. This mechanism is similar to "inertial navigation" or "automatic centering," significantly improving the responsiveness and convenience of cross-section browsing.

[0092] Obtain the corrected rotational angular velocity ω final and its rotation axis direction R axis Then, the system updates the spatial attitude of any current cross-section. The specific operation is as follows: Calculate the single-step rotation angle Δθ=ω final ×Δt, where Δt is the time interval between two adjacent frames (e.g., 1 / 60 of a second).

[0093] The cross-section normal vector n old about the axis of rotation R axis Rotate by an angle Δθ to obtain a new normal vector n. new .

[0094] The in-plane direction vector u old (Cross-section direction) Same rotation axis R axis Rotate by an angle Δθ to obtain a new in-plane direction vector u. new .

[0095] The new normal vector and the in-plane direction vector must remain orthogonal and of unit length. Accuracy can be ensured by re-orthogonalizing (such as Schmidt orthogonalization) or by direct rotation calculation.

[0096] The rotated cross-section is defined as the current initial arbitrary end face output in this step. This cross-section will serve as the reference plane for the sampling ray in the subsequent step S106.

[0097] Application Example: When a user is viewing an abdominal cross-section, the anchor point is the center of the kidney. Initially, the projection of this point onto the cross-sectional image is exactly at the center of the window (ε=0), with no restoring force. The user suddenly rotates the cross-section rapidly to the right, and the kidney center moves to near the right edge of the window (ε≈0.45). The system detects the increased normalized offset distance and generates a restoring angular velocity ω. recovery This is superimposed on the rotational angular velocity input by the user. As the user continues to rotate, they will feel a kind of "resistance" or "returning tendency," and when the user stops rotating (ω... user After the value is 0, the system automatically and slowly rotates the cross-section back to the position near the center of the window in a decaying manner, without requiring the user to manually and precisely correct it.

[0098] In summary, since the angular velocity of the return is proportional to the deviation distance, the return force is stronger when the anchor point deviates far, which can quickly bring the observation focus back to the field of view; when the anchor point is close to the center, the return force weakens, avoiding small oscillations and ensuring the flexibility of fine-tuning the cross section; and the return force only acts on the rotation angle of the cross section and does not change the actual position of the anchor point in three-dimensional space, so the target point of the user positioning will not be lost.

[0099] Step S106: Define a sampling ray for each pixel in the current initial arbitrary end face, and advance point by point along each sampling ray from the starting point. Calculate the absolute difference in gray level between the current sampling point and the previous sampling point and accumulate the difference. Stop sampling the sampling ray when the accumulated value reaches a preset threshold.

[0100] In this step, the direction of the sampling ray is the direction of the normal vector of the current initial arbitrary end face.

[0101] Specifically, for each pixel in the current initial arbitrary end face, the three-dimensional space point corresponding to the pixel is obtained as the starting point of the sampling ray; Obtain the four adjacent pixels of a pixel, including the pixel above, the pixel below, the pixel to the left, and the pixel to the right, and check whether the four adjacent pixels have completed the sampling of the sampling ray and record the sampling stop depth; If at least one neighboring pixel has a recorded sampling stop depth, the arithmetic mean of all recorded sampling stop depths is calculated, and the arithmetic mean is used as the expected stop depth of the pixel; if no neighboring pixels have a recorded sampling stop depth, the preset maximum depth threshold is used as the expected stop depth. Starting from the sampling ray origin, with the normal vector direction as the sampling direction, and using the preset initial step size as the current step size, the process moves forward point by point. During the advancement process, the ratio of the current advanced depth to the expected stopping depth is calculated in real time, and the current step size is dynamically adjusted based on this ratio. When the ratio is less than one-half, keep the current step size as the initial step size; When the ratio is greater than or equal to one-half and less than three-quarters, reduce the current step size to half of the initial step size; When the ratio is greater than or equal to three-quarters, reduce the current step size to one-quarter of the initial step size; At each sampling point, the gray value of the current sampling point is obtained, and the gray value of the previous sampling point is obtained. The absolute difference between the gray value of the current sampling point and the gray value of the previous sampling point is calculated, and the absolute difference is accumulated into the gray value change accumulator. Determine whether the current value of the grayscale change accumulator has reached the preset accumulation threshold: If the current value of the grayscale change accumulator has reached the accumulation threshold, the advancement of the current sampling ray will be stopped immediately, and the current depth advanced will be recorded as the sampling stop depth of the pixel. If the current value of the grayscale change accumulator has not reached the accumulation threshold, it continues to advance with the adjusted step size; If the grayscale change accumulator still has not reached the accumulation threshold when the advancing depth exceeds the preset maximum depth threshold, the advancing will be forcibly stopped, and the maximum depth threshold will be recorded as the sampling stop depth of the pixel.

[0102] In one specific embodiment, for each pixel in the current initial arbitrary end face (the position is determined by the row index i and column index j), the system first obtains the three-dimensional spatial coordinates of the pixel, denoted as the starting point S(i,j). This starting point is located on the cross-sectional plane and is the starting position of the sampling ray in the three-dimensional volume data.

[0103] To accelerate the sampling of the current pixel by utilizing prior information from neighboring pixels that have already been sampled, the system acquires the four neighboring pixels of the current pixel: the pixel above (i-1,j), the pixel below (i+1,j), the pixel to the left (i,j-1), and the pixel to the right (i,j+1). The system sequentially checks whether these four neighboring pixels have completed sampling of the sampling ray and recorded the effective sampling stop depth. Since the sampling order of cross-sectional pixels can follow a certain scanning order (e.g., line-by-line scanning), when processing the current pixel, some neighboring pixels may have already been processed and their depth values ​​recorded in previous iterations. The pixel sampling order can be implemented using a raster scanning order from top to bottom and from left to right, in which case the pixels above and to the left are usually already sampled, while the pixels to the right and to the bottom have not yet been sampled; alternatively, a parallel approach can be used to sample multiple pixels simultaneously, in which case the neighborhood information may depend on the depth estimation of the previous frame.

[0104] If at least one adjacent pixel has a recorded sampling stopping depth, the system calculates the arithmetic mean of all recorded sampling stopping depths and uses this arithmetic mean as the expected stopping depth D of the current pixel. exp This expected value reflects the average position of the neighborhood structure in the depth direction. Utilizing the spatial continuity of the tissue structure in the cross-sectional image, the expected sampling stopping depth of the current pixel should be similar to that of the surrounding pixels.

[0105] If none of the four neighboring pixels of the current pixel have been sampled and their depth recorded (e.g., when processing the first pixel), the system will set the preset maximum depth threshold D. max As the expected stopping depth. D maxThe dimensions of the volume data can be preset, for example, one or two times the diagonal length of the spatial bounding box, to ensure that the sampling ray can penetrate the entire volume data.

[0106] Starting from the sampling ray origin S(i,j), and taking the normal vector direction of the current arbitrary end face (denoted as N) as the sampling direction, the system uses a preset initial step size Δ step0 (For example, 0.5 times the voxel spacing or 0.5 mm) is used as the current step size, and the process moves forward point by point.

[0107] During the advancement process, the system records the current depth d that has been advanced in real time. current (The distance moved along the normal vector direction from the starting point S). Then calculate the ratio of the current propulsion depth to the expected stopping depth: r = d current / D exp (If D) exp If the value is 0, special handling is applied, for example, r=0). Based on the dynamic change of the ratio r, the system adjusts the current step size according to the following rules: When r < 1 / 2, meaning the current depth has not yet reached 50% of the expected depth, and the tissue boundary is still far away, a large stride can be used to advance quickly, maintaining the current stride length as the initial stride length Δ. step0 .

[0108] When 1 / 2 ≤ r < 3 / 4, meaning the current depth has reached 50% but not 75% of the expected depth, and is close to the tissue boundary, the current step size is reduced to half of the initial step size, i.e., Δ. step =Δ step0 / 2, to increase sampling density and prevent missing boundary details.

[0109] When r ≥ 3 / 4, meaning the current depth has reached more than 75% of the expected depth and we have entered the critical region of the tissue boundary, the current step size is further reduced to one-quarter of the initial step size, i.e., Δ. step =Δ step0 / 4, to achieve fine sampling and accurately capture the location of grayscale abrupt changes.

[0110] Step size adjustment is performed before each sampling point advances, ensuring that each step uses a step size that matches the current region. This segmented variable step size strategy is fast and economical when far from the expected boundary, and finely focuses when close to the boundary, ensuring the accuracy of boundary capture while avoiding computational redundancy caused by small step sizes throughout.

[0111] At each sampling point (including every advance position after the starting point), the system performs the following operations: Obtaining grayscale value: Based on the three-dimensional spatial coordinates of the current sampling point, obtain the grayscale value g of that point from the three-dimensional volume data constructed and uploaded to the display memory in step S101 through trilinear interpolation (or nearest neighbor sampling).current .

[0112] Calculate the absolute difference in gray levels: Obtain the gray level value g of the previous sampling point. prev (For the first sampling point, i.e., the gray level at the starting point S, g) prev (Starting with grayscale), calculate the absolute difference Δg = |g current -g prev |

[0113] Accumulation to grayscale change accumulator: Maintain an accumulator variable acc for each pixel, initially set to 0. After each step, accumulate the absolute difference Δg into the accumulator variable acc.

[0114] Stopping condition: The system presets an accumulation threshold T. acc (For example, 50 or 100, depending on the grayscale range and volume data characteristics). Determine if the accumulator variable acc has reached or exceeded T. acc : If the accumulator variable acc ≥ T acc This indicates that the cumulative change in grayscale since sampling began has been large enough, meaning the sampling ray has crossed a clear tissue boundary (e.g., from air into soft tissue, or from soft tissue into bone). At this point, immediately stop the current sampling ray's advance and set the current depth d. current Record the sampling stopping depth D of this pixel. stop(i,j) .

[0115] If the accumulator variable acc <T acc If the adjusted step size is used, then continue to advance one step size forward and repeat the above process.

[0116] Forced Stop: To prevent the ray from never meeting the accumulation threshold (e.g., propagating entirely in a homogeneous medium), the system also sets a maximum depth threshold D. max If the depth d is advanced current Exceeded D max And acc has not yet reached T acc Then, the advancement will be forcibly stopped, and D will be... max This is recorded as the sampling stop depth. This indicates that there is no obvious grayscale abrupt change boundary at this pixel, and the sampling ray penetrates the entire volume data range.

[0117] To maximize the use of neighborhood information, the system can process each pixel in the following priority order: By adopting a raster scanning order (row by row from left to right, from top to bottom), when processing the current pixel, its upper and left neighbors must have already been processed, providing effective depth prior.

[0118] For the boundary pixels of the first row and the first column, some neighborhoods are missing. The system only uses the average depth of the existing neighborhoods. If none of them exist, the maximum depth threshold is used.

[0119] Experiments show that this method can reduce the overall number of sampling points by about 70% to 85%, with minimal impact on boundary detection accuracy.

[0120] Application Example: In a 512×512 abdominal cross-section, a certain pixel corresponds to the abdominal aorta region. The sampling stopping depth of the pixels above and to the left of this pixel has been determined to be 15mm (corresponding to the anterior wall boundary of the aorta). Therefore, the expected stopping depth D is... exp The initial step size is approximately 15mm. The sampling ray originates from the cross-sectional plane (a point inside the aortic lumen) and advances along the normal direction (e.g., from front to back). Initially, the depth is less than 7.5mm (r < 1 / 2), with each step being 0.5mm, and the grayscale accumulator increases slowly. When the depth reaches 7.5mm, the step size is reduced to 0.25mm; when the depth reaches 11.25mm, the step size is reduced to 0.125mm. At a depth of approximately 14.3mm, the grayscale value suddenly changes from low grayscale (blood) in the anterior wall lumen to high grayscale (calcified plaques or adventitia), with a significant increase in single-step Δg. The accumulator instantly reaches the threshold, sampling stops immediately, and the stopping depth of 14.3mm is recorded.

[0121] In summary, by utilizing the spatial continuity of tissue structures in cross-sectional images and using the average depth of adjacent pixels that have been sampled as the expected stopping depth, the search range is significantly narrowed. In areas of abrupt grayscale changes, the step size is automatically switched to a smaller value, which can accurately capture the precise location of tissue boundaries and avoid the problem of missing thin-layer structures or blurry boundaries due to excessively large step sizes.

[0122] Step S107: Resample the volume data resource according to the sampling results of all sampled rays to generate an arbitrary cross-section of the current target, and display the cross-section image corresponding to the arbitrary cross-section of the current target in the arbitrary cross-section view.

[0123] In this step, the sampling stopping depth of all pixels is obtained, and an initial depth map with the same number of rows and columns as the cross-sectional pixel grid is constructed, where the value at each position in the initial depth map corresponds to the sampling stopping depth of the pixel at that position; The initial depth map is subjected to median filtering with a preset window size. The depth value of each pixel is replaced with the median of all depth values ​​within the filtering window to obtain a smooth depth map. For each pixel in the smooth depth map, the depth value of the pixel is obtained, and the rate of change of gray values ​​in the neighborhood of the pixel is obtained. The rate of change is calculated as follows: taking the sampling ray of the pixel as the center, the gray values ​​of adjacent pixels are taken along the cross-sectional plane, and the variance of the gray values ​​taken is used as the rate of change index. The sampling thickness is determined based on the rate of change index, where the sampling thickness is inversely proportional to the rate of change index. The larger the rate of change index, the smaller the sampling thickness, and the smaller the rate of change index, the larger the sampling thickness. The sampling thickness is not less than the preset minimum thickness and not greater than the preset maximum thickness. Centered on the depth value of a pixel, the sampling depth range of the pixel is formed by extending the sampling thickness forward by half and backward by half. Within the sampling depth range, weights are assigned according to the distance of the sampling point from the center depth. The closer the distance, the greater the weight. The weight distribution satisfies the shape of a Gaussian function, where the standard deviation of the Gaussian function is equal to one-third of the sampling thickness. Calculate the weighted average of the gray values ​​of all sampling points within the sampling depth range to obtain the final gray value of the pixel; Arrange the final grayscale values ​​of all pixels according to the row and column order of the cross-sectional pixel grid to generate the initial cross-sectional image; Adaptive contrast stretching of the initial cross-sectional image: Calculate the cumulative distribution function of gray values ​​in the initial cross-sectional image, map the gray value at the first percentile of the cumulative distribution function to zero, map the gray value at the ninety-ninth percentile of the cumulative distribution function to the maximum value, and linearly map the intermediate gray values ​​to obtain the final cross-sectional image; The final cross-sectional image is output to the display window of any cross-sectional view for display.

[0124] In one specific embodiment, the system first reads the sampling stopping depth of each pixel recorded in step S106 (i.e., the advance depth corresponding to the first time the accumulated grayscale change value of the sampling ray reaches a preset threshold, or the maximum depth threshold when forced to stop). These depth values ​​reflect the "effective visible depth" of the tissue structure corresponding to each pixel in the cross-sectional image in the line of sight direction—usually corresponding to the first significant grayscale jump position of the tissue boundary, such as the interface between air and soft tissue, or between soft tissue and bone.

[0125] The system constructs a two-dimensional array with the same number of rows and columns as the cross-sectional pixel grid, called the initial depth map. The value at each position in the array is equal to the sampling stopping depth of the corresponding pixel. Due to the influence of noise or local artifacts during the sampling process, isolated depth outliers may exist in the initial depth map (e.g., the depth of some pixels is significantly different from that of their neighboring pixels). These outliers will introduce noise in subsequent resampling.

[0126] To eliminate outliers, the system performs median filtering on the initial depth map. The filter window size is preset to, for example, 3×3 or 5×5 pixels. For each pixel in the depth map, the depth values ​​of all pixels within the filter window's coverage area are taken, these depth values ​​are sorted by size, and the median (the value in the middle after sorting) is taken as the new depth value for that pixel. Compared to mean filtering, median filtering has better edge-preserving properties, able to remove isolated outliers while retaining true depth boundaries. The filtered depth map is called a smoothed depth map.

[0127] For each pixel in the smoothed depth map, its smoothed depth value represents the approximate location of the tissue boundary where that pixel is situated. However, taking only the grayscale value at a single depth location as the final pixel value can easily lead to information loss (e.g., thin layer structures or partial volume effects). Therefore, weighted fusion is needed within a range around that depth to improve image quality and noise resistance. Since the complexity of tissue structures varies in different regions, the fusion range (sampling thickness) should be adaptively adjusted: in regions with rich texture and dramatic grayscale changes (such as bone trabeculae and organ boundaries), a smaller sampling thickness should be used to preserve details; in flat and uniform regions (such as muscle parenchyma and blood pools), a larger sampling thickness can be used to suppress noise.

[0128] The system takes the sampling ray of the current pixel as the center and samples the gray values ​​of neighboring pixels along the cross-sectional plane (i.e., the row and column directions of the cross-sectional image). It calculates the variance of these gray values ​​as a rate of change indicator. A larger variance indicates more drastic local gray-level changes and a more complex organizational structure; a smaller variance indicates a flatter and more uniform local area.

[0129] Then, the system determines the sampling thickness of each pixel based on the rate of change index. The sampling thickness is inversely proportional to the rate of change index: the larger the rate of change index, the smaller the sampling thickness; the smaller the rate of change index, the larger the sampling thickness. Specifically, a linear mapping can be used: let the rate of change index be σ², and the preset minimum thickness be T. min (e.g., 0.5 mm), maximum thickness is T max (For example, 3 mm), then the sampling thickness T = T max -(T max -T min )×(σ²-σ² min ) / (σ² max -σ² min ), where σ² min and σ² max These are the minimum and maximum values ​​of the rate of change index across the entire cross-section, respectively. Meanwhile, T is restricted to [T...]. min ,T max Within the range.

[0130] For each pixel, the depth value of that point in the smoothed depth map is taken as the center depth D. center Using the sampling thickness T calculated in the previous step as the interval radius, extend it forward by T / 2 and backward by T / 2 to form the sampling depth interval [D] of that pixel. center -T / 2, D center +T / 2]. This interval covers a thin layer centered on the organizational boundary, containing grayscale information of both internal and external organizations.

[0131] Within this sampling depth range, the system needs to assign weights based on the distance of each sampling point from the center depth to calculate a weighted average. The weight distribution adopts a Gaussian function shape, meaning that sampling points closer to the center depth have larger weights, and those farther away have smaller weights. Specific parameters include the standard deviation σ of the Gaussian function. gauss Setting it to T / 3 makes the weight at the interval endpoints approximately 5% of the peak value, which can effectively suppress far-field noise.

[0132] Suppose that the sampling ray contains N sampling points within the depth interval (obtained from the advancement result of step S106), and the depth of each sampling point k is depth. k The grayscale value is g k Then the weight w k =exp(-(depth k -D center ) 2 / (2×(T / 3) 2 These weights are normalized (i.e., divided by the sum of the weights of all sampling points), and then the weighted average gray value is calculated.

[0133] The weighted average value is used as the final grayscale value of the pixel. The advantage of using Gaussian weighting is that it can fully preserve the dominant information of the tissue boundary at the center depth, while gently blending the grayscale of the neighboring areas, reducing image flicker caused by small fluctuations in sampling depth, and avoiding the blurring of boundaries caused by simple averaging.

[0134] Arrange the final grayscale values ​​of all pixels according to the row and column order of the cross-sectional pixel grid (i.e., row index from top to bottom, column index from left to right) to generate a two-dimensional grayscale matrix, called the initial cross-sectional image. Since the original grayscale range of volume data (such as HU units for CT values) is often very wide (-1024~3071), when directly mapped to a display device, the effective information may be concentrated in a certain sub-interval, resulting in low image contrast. For example, when observing bones, the grayscale values ​​of soft tissue and air occupy a large portion of the dynamic range, making the bones themselves appear darker.

[0135] Therefore, the system performs adaptive contrast stretching on the initial cross-sectional image. The specific method is as follows: The grayscale value distribution of all pixels in the initial cross-sectional image is statistically analyzed, and the cumulative distribution function (CDF) of this distribution is calculated. The CDF is a monotonically increasing function, and F(g) represents the proportion of pixels with grayscale values ​​≤ g.

[0136] Set the low truncation percentile (e.g., 1%) and the high truncation percentile (e.g., 99%). Find the grayscale value g. low , so that F(g) low =1%; Search for grayscale value g high , so that F(g) high )=99%. Thus, the grayscale value is lower than g. low Pixels with gray values ​​higher than g will be considered abnormal dark spots (such as noise or air). high Pixels that are considered abnormal bright spots (such as metallic artifacts) will be truncated to expand the contrast of the main grayscale range.

[0137] g low Mapped to the minimum brightness value of the display device (usually 0), g high Mapped to the maximum brightness value (e.g., 255 or 65535). For grayscale values ​​g between g low and g high The pixels between them are linearly mapped: g mapped =(gg low ) / (g high -g low )×(max display -min display )+min display , where g mapped The maximum value is the final brightness value of the pixel output to the display window after contrast adaptive stretching. display min is the maximum brightness value that the display window can output. display The maximum brightness value that the display window can output; After all pixels are mapped as described above, the final cross-sectional image is obtained.

[0138] Contrast adaptive stretching can effectively utilize the dynamic range of the display device, making the tissue structure in cross-sectional images more vivid and the details clearer.

[0139] The final cross-sectional image is output as a texture map or direct bitmap to the display window of any cross-sectional view for real-time display. Since the entire calculation process (median filtering, sampling thickness calculation, Gaussian weighting, contrast stretching) can be efficiently implemented in parallel on a graphics processing unit (GPU) or a central processing unit (CPU), the refresh rate can reach tens of frames per second, supporting continuous user interaction.

[0140] Application Example: Taking an arbitrary cross-section of an abdominal CT scan as an example, the cross-sectional image covers structures such as the lumbar spine, abdominal aorta, and kidneys. The original sampling stopping depth map may contain isolated noise (such as a pixel stopping prematurely due to intestinal gas interference, resulting in a shallow depth). Median filtering can effectively eliminate such isolated points. In the bone-soft tissue boundary region, the variance of the grayscale change rate is large, so the sampling thickness is adaptively reduced to 0.8mm. After Gaussian weighting, the edge of the bone cortex in this region is clear and sharp. In the flat area inside the liver parenchyma, the variance is small, so the sampling thickness is automatically increased to 2.5mm. Gaussian weighting makes the grayscale of the liver tissue smoother and more uniform, reducing speckle noise. Finally, after contrast stretching, the originally dark lumbar spine and kidney structures become prominent, and the grayscale difference between abdominal fat and muscle is amplified, allowing doctors to easily identify anatomical details.

[0141] Please see Figure 2 The diagram illustrates the structural block diagram of a target-point driven arbitrary cross-section CT image visualization system according to this application.

[0142] like Figure 2 As shown, the arbitrary cross-section CT image visualization system 200 includes an acquisition module 210, a mapping module 220, a calculation module 230, a generation module 240, an update module 250, a sampling module 260, and an output module 270.

[0143] The acquisition module 210 is configured to acquire a target medical image sequence, construct three-dimensional volume data based on the spatial location information and grayscale values ​​of each frame of the target medical image, and upload the three-dimensional volume data to the graphics processor's video memory to form a volume data resource that can be shared by any cross-sectional view; the mapping module 220 is configured to receive the selected point coordinates input by the user and map the selected point coordinates to the three-dimensional volume data according to a preset coordinate mapping strategy to obtain an initial target point; the calculation module 230 is configured to take the initial target point as the starting mass, use the grayscale value of each voxel in the three-dimensional volume data as a virtual mass, calculate the virtual gravitational resultant force of all voxels in a preset spherical neighborhood centered on the current position on the starting mass, and move the starting mass along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than a preset gravitational threshold, and define the final stopping position as the anchoring target point; the generation module 240 is configured to take the anchoring target point as the origin of the tangent plane, and calculate the virtual gravitational resultant force of all voxels in a preset spherical neighborhood centered on the current position on the starting mass, and move the starting mass along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than a preset gravitational threshold, and define the final stopping position as the anchoring target point; The module 250 generates an arbitrary cross-section based on the initial normal vector and the in-plane direction vector; the update module 250 is configured to obtain the rotational angular velocity input by the user, generate an inertia coefficient based on the offset vector between the projection point of the anchored target point on the current arbitrary cross-section and the geometric center of the current arbitrary cross-section, correct the rotational angular velocity using the inertia coefficient, and update the spatial attitude of the current arbitrary cross-section to obtain the current initial arbitrary end face; the sampling module 260 is configured to define a sampling ray for each pixel in the current initial arbitrary end face, advance along each sampling ray from the starting point, calculate the absolute difference in grayscale between the current sampling point and the previous sampling point in turn, and accumulate the difference. When the accumulated value reaches a preset threshold, the sampling of the sampling ray is stopped; the output module 270 is configured to resample the volume data resource based on the sampling results of all sampling rays, generate the current target arbitrary cross-section, and display the cross-section image corresponding to the current target arbitrary cross-section in the arbitrary cross-section view.

[0144] It should be understood that Figure 2 The modules and references described in the document Figure 1 The steps described in the text correspond to those in the method described above. Therefore, the operations, features, and corresponding technical effects described above also apply to the method described in the text. Figure 2 The various modules in the document will not be described in detail here.

[0145] In other embodiments, the present invention also provides a computer-readable storage medium having a computer program stored thereon, wherein when the program instructions are executed by a processor, the processor performs the target point-driven arbitrary section CT image visualization method in any of the above method embodiments. In one embodiment, the computer-readable storage medium of the present invention stores computer-executable instructions, which are configured as follows: The target medical image sequence is acquired, and three-dimensional volume data is constructed based on the spatial location information and grayscale value of each frame of the target medical image. The three-dimensional volume data is then uploaded to the video memory of the graphics processor to form a volume data resource that can be shared by any cross-sectional view. Receive the selected point coordinates input by the user, and map the selected point coordinates to the three-dimensional volume data according to the preset coordinate mapping strategy to obtain the initial target point; Using the initial target point as the starting mass, the gray value of each voxel in the three-dimensional volume data is used as the virtual mass. The virtual gravitational resultant force of all voxels in the preset spherical neighborhood centered on the current position is calculated. The starting mass is moved along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold. The final stopping position is defined as the anchoring target point. Using the anchored target point as the origin of the tangent, generate the current arbitrary cross section based on the initial normal vector and the in-plane direction vector; The rotational angular velocity input by the user is obtained. An inertia coefficient is generated based on the offset vector between the projection point of the anchored target point on the current arbitrary cross-section and the geometric center of the current arbitrary cross-section. The rotational angular velocity is corrected using the inertia coefficient, and the spatial attitude of the current arbitrary cross-section is updated to obtain the current initial arbitrary end face. Define a sampling ray for each pixel in the current initial arbitrary end face, and advance point by point along each sampling ray from the starting point. Calculate the absolute difference in gray level between the current sampling point and the previous sampling point and accumulate the difference. Stop sampling the sampling ray when the accumulated value reaches a preset threshold. The volume data resource is resampled based on the sampling results of all sampled rays to generate an arbitrary cross-section of the current target, and the cross-sectional image corresponding to the arbitrary cross-section of the current target is displayed in the arbitrary cross-sectional view.

[0146] Computer-readable storage media may include a program storage area and a data storage area, wherein the program storage area may store an operating system and an application program required for at least one function; the data storage area may store data created based on the use of the target-point driven arbitrary-section CT image visualization system, etc. Furthermore, the computer-readable storage medium may include high-speed random access memory, and may also include memory, such as at least one disk storage device, flash memory device, or other non-volatile solid-state storage device. In some embodiments, the computer-readable storage medium may optionally include memory remotely disposed relative to a processor, which can be connected to the target-point driven arbitrary-section CT image visualization system via a network. Examples of such networks include, but are not limited to, the Internet, corporate intranets, local area networks, mobile communication networks, and combinations thereof.

[0147] Figure 3This is a schematic diagram of the structure of the electronic device provided in the embodiment of the present invention, such as... Figure 3 As shown, the device includes a processor 310 and a memory 320. The electronic device may also include an input device 330 and an output device 340. The processor 310, memory 320, input device 330, and output device 340 can be connected via a bus or other means. Figure 3 Taking a bus connection as an example, the memory 320 is the computer-readable storage medium described above. The processor 310 executes various server functions and data processing by running non-volatile software programs, instructions, and modules stored in the memory 320, thereby implementing the target-point-driven arbitrary section CT image visualization method described in the above embodiment. The input device 330 can receive input digital or character information and generate key signal inputs related to user settings and function control of the target-point-driven arbitrary section CT image visualization system. The output device 340 may include a display screen or other display device.

[0148] The aforementioned electronic device can execute the method provided in the embodiments of the present invention, and has the corresponding functional modules and beneficial effects for executing the method. Technical details not described in detail in this embodiment can be found in the method provided in the embodiments of the present invention.

[0149] In one implementation, the above-described electronic device is applied to a target-point driven arbitrary section CT image visualization system for a client, comprising: at least one processor; and a memory communicatively connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to enable the at least one processor to: The target medical image sequence is acquired, and three-dimensional volume data is constructed based on the spatial location information and grayscale value of each frame of the target medical image. The three-dimensional volume data is then uploaded to the video memory of the graphics processor to form a volume data resource that can be shared by any cross-sectional view. Receive the selected point coordinates input by the user, and map the selected point coordinates to the three-dimensional volume data according to the preset coordinate mapping strategy to obtain the initial target point; Using the initial target point as the starting mass, the gray value of each voxel in the three-dimensional volume data is used as the virtual mass. The virtual gravitational resultant force of all voxels in the preset spherical neighborhood centered on the current position is calculated. The starting mass is moved along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold. The final stopping position is defined as the anchoring target point. Using the anchored target point as the origin of the tangent, generate the current arbitrary cross section based on the initial normal vector and the in-plane direction vector; The rotational angular velocity input by the user is obtained. An inertia coefficient is generated based on the offset vector between the projection point of the anchored target point on the current arbitrary cross-section and the geometric center of the current arbitrary cross-section. The rotational angular velocity is corrected using the inertia coefficient, and the spatial attitude of the current arbitrary cross-section is updated to obtain the current initial arbitrary end face. Define a sampling ray for each pixel in the current initial arbitrary end face, and advance point by point along each sampling ray from the starting point. Calculate the absolute difference in gray level between the current sampling point and the previous sampling point and accumulate the difference. Stop sampling the sampling ray when the accumulated value reaches a preset threshold. The volume data resource is resampled based on the sampling results of all sampled rays to generate an arbitrary cross-section of the current target, and the cross-sectional image corresponding to the arbitrary cross-section of the current target is displayed in the arbitrary cross-sectional view.

[0150] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., including several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods of various embodiments or some parts of embodiments.

[0151] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A target-point driven method for visualizing arbitrary cross-sectional CT images, characterized in that, include: The target medical image sequence is acquired, and three-dimensional volume data is constructed based on the spatial location information and grayscale value of each frame of the target medical image. The three-dimensional volume data is then uploaded to the video memory of the graphics processor to form a volume data resource that can be shared by any cross-sectional view. Receiving user-inputted coordinates of selected points and mapping them to the 3D volume data according to a preset coordinate mapping strategy to obtain an initial target point, including: The input time of the selected point coordinates is obtained, and the coordinates of multiple consecutive historical selected points before the input time are recorded to form a coordinate sequence; The user's line of sight is fitted on the display window plane according to the coordinate sequence, and the instantaneous velocity direction and instantaneous acceleration direction of the movement trajectory at the current moment are calculated, wherein the instantaneous velocity direction and instantaneous acceleration direction are both two-dimensional direction vectors in the display window plane; Starting from the selected point coordinates, a ray is projected into the display window plane along the instantaneous velocity direction, and the ray is extended from the display window into the three-dimensional scene to form a first spatial ray. The first intersection point of the first spatial ray and the spatial bounding box of the three-dimensional volume data is calculated as the first intersection point, and the distance between the first intersection point and the observation point is used as the first depth candidate value. The observation point is the location of the virtual camera in the three-dimensional scene. Starting from the selected point coordinates, another ray is projected into the display window plane along the instantaneous acceleration direction, and the other ray is extended from the display window into the three-dimensional scene to form a second spatial ray. The first intersection point of the second spatial ray and the spatial bounding box is calculated as the second intersection point, and the distance between the second intersection point and the observation point is used as the second depth candidate value. Calculate the angle between the instantaneous velocity direction and the instantaneous acceleration direction, and determine whether the angle is less than a preset angle threshold; If the included angle is less than a preset angle threshold, the weighted average of the first depth candidate value and the second depth candidate value is used as the initial depth value of the initial target point. If the included angle is not less than a preset angle threshold, then the first depth candidate value is used as the initial depth value of the initial target point; The initial target point is constructed by using the selected point coordinates as two-dimensional plane coordinates and the initial depth value as spatial coordinates; Using the initial target point as the starting mass, the gray value of each voxel in the three-dimensional volume data is used as the virtual mass. The virtual gravitational resultant force of all voxels in the preset spherical neighborhood centered on the current position is calculated. The starting mass is moved along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold. The final stopping position is defined as the anchoring target point. Using the anchored target point as the origin of the tangent, generate the current arbitrary cross section based on the initial normal vector and the in-plane direction vector; The rotational angular velocity input by the user is obtained. An inertia coefficient is generated based on the offset vector between the projection point of the anchored target point on the current arbitrary cross-section and the geometric center of the current arbitrary cross-section. The rotational angular velocity is corrected using the inertia coefficient, and the spatial attitude of the current arbitrary cross-section is updated to obtain the current initial arbitrary cross-section. Define a sampling ray for each pixel in the current initial arbitrary cross section, and advance point by point along each sampling ray from the starting point. Calculate the absolute difference in gray level between the current sampling point and the previous sampling point and accumulate the difference. Stop sampling the sampling ray when the accumulated value reaches a preset threshold. The volume data resource is resampled based on the sampling results of all sampled rays to generate an arbitrary cross-section of the current target, and the cross-sectional image corresponding to the arbitrary cross-section of the current target is displayed in the arbitrary cross-sectional view.

2. The target-point driven arbitrary section CT image visualization method according to claim 1, characterized in that, The step of taking the initial target point as the starting mass point, using the gray value of each voxel in the three-dimensional volume data as the virtual mass, and calculating the virtual gravitational resultant force of all voxels in a preset spherical neighborhood centered on the current position on the starting mass point includes: The particle at the current iteration position is taken as the current particle, wherein the current particle is the starting particle during the first iteration; A spherical neighborhood is constructed with the current position of the current particle as the center point, and the radius of the spherical neighborhood is a preset initial radius; Obtain the gray values ​​of all voxels within the spherical neighborhood, and statistically analyze the distribution of the gray values ​​to determine the first gray value with the highest frequency and the second gray value with the second highest frequency. Calculate the absolute difference between the first gray value and the second gray value, and determine whether the absolute difference is greater than a preset separation threshold; If the absolute difference is greater than the preset separation threshold, it is determined that there are at least two different tissue types in the spherical neighborhood, and the gray centroid corresponding to each tissue type is calculated respectively. Voxels with gray values ​​within the first preset gray value range are classified as the first tissue type, and voxels with gray values ​​within the second preset gray value range are classified as the second tissue type. Calculate the first gravitational vector pointing from the current particle to the gray-scale centroid of the first tissue type, and the second gravitational vector pointing from the current particle to the gray-scale centroid of the second tissue type, respectively. The magnitude of the first gravitational vector is proportional to the proportion of the number of voxels contained in the first tissue type to the total number of voxels in the spherical neighborhood, and the magnitude of the second gravitational vector is proportional to the proportion of the number of voxels contained in the second tissue type to the total number of voxels in the spherical neighborhood. The first gravitational vector and the second gravitational vector are vector superimposed to obtain the virtual gravitational resultant force on the current particle. If the absolute difference is not greater than a preset separation threshold, the spherical neighborhood is determined to be a single tissue type. The gray value of each voxel in the spherical neighborhood is normalized to the range of zero to one and used as the virtual mass of the voxel. The gravitational vector of each voxel to the center point is then calculated. All gravitational vectors are summed to obtain the virtual gravitational resultant force. The direction of the gravitational vector of each voxel is the direction from the voxel to the center point. The magnitude of the gravitational vector is directly proportional to the virtual mass of the voxel and inversely proportional to the square of the distance from the voxel to the center point.

3. The target-point driven arbitrary section CT image visualization method according to claim 1, characterized in that, Moving the starting particle along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than a preset gravitational threshold, and defining the final stopping position as the anchoring target point, includes: For the current particle at the current iteration position, where the current particle is the starting particle during the first iteration, obtain the current position of the current particle and the current virtual gravitational resultant force at the current position; Determine whether the magnitude of the current virtual gravitational resultant force is less than a preset gravitational threshold; If the force is less than the preset gravity threshold, the movement stops immediately, and the position of the current particle is defined as the anchoring target point. If it is not less than the preset gravity threshold, then obtain the historical position before the last movement, and calculate the historical displacement vector from the historical position to the current mass point; Calculate the angle between the historical displacement vector and the current virtual gravitational resultant force direction, and determine whether the angle is greater than a preset oscillation angle threshold in a series of preset iterations; If the included angle is greater than the preset oscillation angle threshold in consecutive preset iterations, it is determined that the current particle is trapped in a local oscillation state, and an oscillation suppression operation is performed to obtain a new particle position. Otherwise, the direction of the current virtual gravitational resultant force is directly used as the moving direction, and the current particle is moved with the current moving step size to obtain a new particle position. The current moving step size is proportional to the magnitude of the current virtual gravitational resultant force. The new particle position is used as the current particle position in the next iteration. The virtual gravitational resultant force is recalculated until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold. At this point, the position of the current particle is defined as the anchoring target point.

4. The target-point driven arbitrary section CT image visualization method according to claim 3, characterized in that, The process of performing oscillation suppression to obtain new particle positions includes: Record the oscillation axis direction, which is the average direction of multiple consecutive historical displacement vectors; The current virtual gravitational force is decomposed into a first component along the oscillation axis and a second component perpendicular to the oscillation axis. Multiply the first component by a preset damping coefficient to obtain the suppressed first component; The suppressed first component and the second component are recombined to obtain the suppressed virtual gravitational resultant force; Using the direction of the suppressed virtual gravitational resultant force as the direction of movement, and taking half of the current movement step as the new movement step, the current particle is moved along the direction of movement to obtain a new particle position.

5. The target-point driven arbitrary section CT image visualization method according to claim 1, characterized in that, in, The current arbitrary cross-section contains a two-dimensional pixel grid determined by pixel width, pixel height and sampling interval, and each pixel in the two-dimensional pixel grid corresponds to a three-dimensional spatial point in the three-dimensional volume data; The step of generating any current cross-section based on the initial normal vector and the in-plane direction vector, with the anchored target point as the origin of the tangent, includes: Obtain the local neighborhood of the anchor target point, wherein the local neighborhood is a cubic region with a preset side length centered at the anchor target point; Calculate the gray-level gradient of each voxel in the local neighborhood to obtain the gradient direction vector of each voxel, and count the direction distribution of all gradient direction vectors to construct a gradient direction histogram. The horizontal axis of the gradient direction histogram is the direction angle interval, and the vertical axis is the number of gradient direction vectors falling within the direction angle interval. The principal gradient direction and the secondary gradient direction are identified from the gradient direction histogram, wherein the principal gradient direction is the direction corresponding to the direction angle interval with the highest frequency of occurrence, and the secondary gradient direction is the direction with the highest frequency of occurrence in the direction angle interval orthogonal to the principal gradient direction; The subgradient direction is used as a candidate direction for the cross-sectional normal vector, and the principal gradient direction is used as a candidate direction for the in-plane direction vector. Calculate the target angle between the candidate direction and the initial normal vector input by the user; If the target angle is less than a preset first angle threshold, then the candidate direction is used as the normal vector of any current cross section. If the target included angle is greater than a preset second angle threshold, then the initial normal vector is used as the normal vector of any current cross section, wherein the second angle threshold is greater than the first angle threshold; If the included angle is between the first angle threshold and the second angle threshold, then the weighted average direction of the candidate direction and the initial normal vector is calculated as the normal vector of the current arbitrary section. Determine the target direction that is simultaneously perpendicular to the normal vector of the current arbitrary cross section and the preset global upward direction vector. Take the target direction as the in-plane direction vector and normalize the in-plane direction vector to obtain the normalized in-plane direction vector. The global upward direction vector is a unit vector pointing vertically upward in three-dimensional space. Using the anchored target point as the origin of the tangent, and based on the normal vector of the current arbitrary cross-section and the normalized in-plane direction vector, combined with the size of the two-dimensional pixel grid and the sampling interval, the coordinates of the three-dimensional spatial point corresponding to each pixel are calculated to generate the current arbitrary cross-section.

6. The target-point driven arbitrary section CT image visualization method according to claim 1, characterized in that, The step of generating an inertia coefficient based on the offset vector between the projection point of the anchored target point on the current arbitrary cross-section and the geometric center of the current arbitrary cross-section, using the inertia coefficient to correct the rotational angular velocity, and updating the spatial attitude of the current arbitrary cross-section to obtain the current initial arbitrary cross-section includes: Calculate the vertical projection of the anchoring target point onto the plane of the current arbitrary section; Obtain the geometric center point of the current arbitrary cross-section display window, calculate the offset vector from the geometric center point to the vertical projection point, and calculate the magnitude of the offset vector; Obtain the diagonal length of the display window and calculate the normalized offset distance; A one-dimensional virtual spring-mass system is constructed with the geometric center point as the origin and the direction pointed to by the offset vector as the positive direction. The vertical projection point is regarded as a mass connected to the spring, the geometric center point is the fixed point at the other end of the spring, and the spring constant is proportional to the normalized offset distance. According to Hooke's Law, the restoring force on the particle is calculated, and the restoring acceleration is calculated based on the preset virtual particle mass. Obtain the scalar value of the rotational angular velocity input by the user, and obtain the rotational axis direction input by the user, wherein the rotational axis direction is perpendicular to the current arbitrary cross-section and passes through the geometric center of the display window; The recovery acceleration is converted into recovery angular velocity in the following manner: the magnitude of the recovery angular velocity is equal to the recovery acceleration divided by the diagonal length of the display window, and the direction of the recovery angular velocity is consistent with the direction of the rotation axis. The corrected rotational angular velocity is obtained by weighted fusion of the recovery angular velocity and the rotational angular velocity input by the user. Based on the corrected rotational angular velocity and the rotation axis direction, the rotation update amount of the normal vector and in-plane direction vector of the current arbitrary cross-section is calculated to obtain the new spatial attitude, and the updated cross-section is defined as the current initial arbitrary cross-section.

7. A target-point driven arbitrary section CT image visualization system, characterized in that, include: The acquisition module is configured to acquire the target medical image sequence, construct three-dimensional volume data based on the spatial location information and grayscale value of each frame of the target medical image, and upload the three-dimensional volume data to the graphics processor's video memory to form a volume data resource that can be shared by any cross-sectional view. The mapping module is configured to receive user-inputted coordinates of selected points and map these coordinates to the 3D volume data according to a preset coordinate mapping strategy to obtain an initial target point, including: The input time of the selected point coordinates is obtained, and the coordinates of multiple consecutive historical selected points before the input time are recorded to form a coordinate sequence; The user's line of sight is fitted on the display window plane according to the coordinate sequence, and the instantaneous velocity direction and instantaneous acceleration direction of the movement trajectory at the current moment are calculated, wherein the instantaneous velocity direction and instantaneous acceleration direction are both two-dimensional direction vectors in the display window plane; Starting from the selected point coordinates, a ray is projected into the display window plane along the instantaneous velocity direction, and the ray is extended from the display window into the three-dimensional scene to form a first spatial ray. The first intersection point of the first spatial ray and the spatial bounding box of the three-dimensional volume data is calculated as the first intersection point, and the distance between the first intersection point and the observation point is used as the first depth candidate value. The observation point is the location of the virtual camera in the three-dimensional scene. Starting from the selected point coordinates, another ray is projected into the display window plane along the instantaneous acceleration direction, and the other ray is extended from the display window into the three-dimensional scene to form a second spatial ray. The first intersection point of the second spatial ray and the spatial bounding box is calculated as the second intersection point, and the distance between the second intersection point and the observation point is used as the second depth candidate value. Calculate the angle between the instantaneous velocity direction and the instantaneous acceleration direction, and determine whether the angle is less than a preset angle threshold; If the included angle is less than a preset angle threshold, the weighted average of the first depth candidate value and the second depth candidate value is used as the initial depth value of the initial target point. If the included angle is not less than a preset angle threshold, then the first depth candidate value is used as the initial depth value of the initial target point; The initial target point is constructed by using the selected point coordinates as two-dimensional plane coordinates and the initial depth value as spatial coordinates; The calculation module is configured to take the initial target point as the starting mass, use the gray value of each voxel in the three-dimensional volume data as the virtual mass, calculate the virtual gravitational resultant force of all voxels in the preset spherical neighborhood centered on the current position on the starting mass, and move the starting mass along the direction of the virtual gravitational resultant force until the magnitude of the virtual gravitational resultant force is less than the preset gravitational threshold, and define the final stopping position as the anchoring target point. The generation module is configured to use the anchored target point as the origin of the tangent and generate any current cross section based on the initial normal vector and the in-plane direction vector. The update module is configured to obtain the rotational angular velocity input by the user, generate an inertia coefficient based on the offset vector between the projection point of the anchored target point on the current arbitrary cross-section and the geometric center of the current arbitrary cross-section, use the inertia coefficient to correct the rotational angular velocity, and update the spatial attitude of the current arbitrary cross-section to obtain the current initial arbitrary cross-section. The sampling module is configured to define a sampling ray for each pixel in the current initial arbitrary cross section, advance along each sampling ray from the starting point point, calculate the absolute difference in gray level between the current sampling point and the previous sampling point in turn and accumulate it, and stop sampling the sampling ray when the accumulated value reaches a preset threshold. The output module is configured to resample the volume data resource based on the sampling results of all sampled rays, generate an arbitrary cross-section of the current target, and display the cross-sectional image corresponding to the arbitrary cross-section of the current target in the arbitrary cross-sectional view.

8. An electronic device, characterized in that, include: At least one processor, and a memory communicatively connected to the at least one processor, wherein the memory stores instructions executable by the at least one processor to enable the at least one processor to perform the method according to any one of claims 1 to 6.

9. A computer-readable storage medium having a computer program stored thereon, characterized in that, When the program is executed by a processor, it implements the method described in any one of claims 1 to 6.

Images