CN101422352B - Interactive coronary artery virtual angioscope implementation method - Google Patents

Abstract

The invention discloses a method for realizing an interactive coronary artery virtual angioscope, which belongs to the technical field of medical imaging. The purpose is to solve the problem of visualization diagnosis and treatment of coronary artery. The technical solution is: it obtains the three-dimensional model of the blood vessel by combining the three-dimensional geometric shape information of the lumen obtained from the approximately orthogonal X-ray coronary angiography image with the cross-sectional data of the lumen obtained by intravascular ultrasound, and then uses virtual The Reality Modeling Language interactively describes the vascular model, enabling visualization of coronary arteries in endoscopic walk-through mode. The invention realizes the interactive visit and display of the three-dimensional blood vessel model, and provides an ideal solution for the development of coronary atherosclerotic lesions, the visual diagnosis and treatment of coronary heart disease, the evaluation of interventional treatment effects, and the training of medical personnel. platform.

Description

A Realization Method of Interactive Coronary Artery Virtual Angioscope

technical field

The invention relates to a realization method of an interactive coronary virtual angioscope based on fusion of multiple imaging methods, which belongs to the technical field of medical imaging.

Background technique

Virtual Endoscope technology (Virtual Endoscope, VE) uses medical images as raw data, comprehensively utilizes digital image processing, computer graphics, scientific computing visualization, virtual reality and other technologies to reconstruct three-dimensional images to form virtual human tissues; Put it into the reconstructed organ cavity, and use navigation or roaming technology and pseudo-color technology to realistically simulate cavity endoscopy.

Coronary Artery Angiography (CAG, Coronary Artery Angiography) and Intravascular Ultrasound (IVUS, Intravascular Ultrasound) are commonly used clinically for coronary artery intervention imaging. The position of the ultrasound probe in the lumen and the structure of the corresponding vessel wall. CAG and IVUS have the characteristics of complementary advantages and disadvantages: CAG reflects the projection contour of the vascular cavity filled with contrast agent, can diagnose diseases such as ischemic heart disease and coronary artery malformation, and is effective for intracoronary thrombolysis, PTCA (percutaneous Interventional surgery such as intraluminal coronary angioplasty) is of great significance, but it cannot provide structural information of the vessel wall and the degree of lesion; IVUS can clearly display the cross-section of the vessel, perform histological classification according to the acoustic characteristics of the plaque, and find the plaques that cannot be displayed by CAG Vascular lesions, observe bifurcations or blurry lesions where vessels overlap, etc. However, due to the use of high-frequency ultrasonic probes, the detection depth is affected, and only a certain section of diseased blood vessels can be measured, and the severely narrowed lumen cannot be entered, and the axial position and spatial direction of the cross-section cannot be determined.

In addition, interventional imaging detection also includes coronary angioscopy, which is a tiny endoscopic imaging technology using fiber optic technology. However, this technique has not been widely accepted clinically. The reasons include: it can only provide morphological data on the surface of the lumen, and cannot observe the deep structure of the lesion in the vessel wall, nor can it perform quantitative analysis of the degree of stenosis and blood flow; It is used to visualize the lesion at the opening of the aorta-coronary artery and the proximal end of the anterior descending and circumflex branches; the blood flow entering from the side hole will blur the field of vision; the lack of maneuverability of the catheter limits the scope of imaging; check The process needs to temporarily block the blood flow, which will lead to the occurrence of myocardial ischemia and so on.

At present, non-invasive cardiovascular imaging mainly includes CTA (CT Angiography) and MRCA (Magnetic Resonance Coronary Angiography). However, the main limitation of cardiac CT examination is that it is easy to produce artifacts and affect the image quality. For MRCA, because the coronary vessels themselves are thin, tortuous and complex in structure, and have the influence of heartbeat and respiration, signals such as fat tissue and myocardial tissue around the coronary arteries can affect the imaging results. At the same time, although there is no radiation during the MR examination process, it is relatively safe, but it is affected by noise, and the safety of some metal implants (such as artificial metal valves, cardiac pacemakers, etc.) has also attracted attention. In short, due to the deficiency and technical limitations caused by the imaging principle, so far CTA and MRCA can only be used as a selection method for the comprehensive evaluation of the heart, or as a screening measure for coronary heart disease catheter angiography to reduce unnecessary Necessary invasive examinations, as well as non-invasive follow-up studies on the effects of cardiac surgery or interventional therapy, cannot completely replace interventional imaging methods in the clinical diagnosis and treatment of coronary heart disease.

In summary, interventional CAG and IVUS are still the main imaging methods for clinical diagnosis and treatment of coronary heart disease, and the two have complementary advantages and disadvantages. At present, there is no interactive coronary virtual angioscopy system based on CAG and IVUS image fusion, which can realize coronary artery visualization in endoscopic roaming mode.

Contents of the invention

The purpose of the present invention is to overcome the deficiencies of the prior art and provide a method for realizing an interactive coronary virtual angioscope capable of realizing coronary artery visualization in endoscopic roaming mode.

The said problem of the present invention is realized with following technical scheme:

A realization method of an interactive coronary virtual angioscope, which fuses the three-dimensional geometric shape information of the lumen obtained from an approximately orthogonal X-ray coronary angiography image with the cross-sectional data of the lumen obtained by intravascular ultrasound, Obtain the three-dimensional model of the blood vessel, and then use the virtual reality modeling language (VRML) to interactively describe the blood vessel model to realize the visualization of the coronary artery in the endoscopic roaming mode. The specific steps are as follows:

A. Simultaneous acquisition of intravascular ultrasound IVUS and X-ray coronary angiography CAG images of the vessel segment of interest:

Place the mechanical ultrasound catheter probe at the distal end of the vessel segment of interest, and use the intravascular ultrasound imager to acquire isometric IVUS image sequences in the manner of ECG gating during the uniform and equidistant retraction of the guide wire , that is, the R wave of the ECG signal is used as a trigger, and images are only collected at the same cardiac phase in each cardiac cycle, which can solve the problem of motion artifacts in coronary IVUS image sequences. At the same time, use the C-arm single-sided X-ray angiography machine to shoot two approximately vertical CAG images at the starting point of the catheter retraction path to record the same heart state;

B, using the collected IVUS and CAG images to establish a three-dimensional model of blood vessels;

C. The virtual reality modeling language is used to realize the visualization of the endoscopic walk-through mode of the coronary artery reconstruction results.

The implementation method of the above-mentioned interactive coronary artery virtual angioscope, the specific steps of using the collected IVUS and CAG images to establish a three-dimensional model of blood vessels are as follows:

a. Based on two CAG images at approximately vertical angles, three-dimensionally reconstruct the retraction path of the ultrasound catheter;

b. Three-dimensional reconstruction of the vessel lumen from the CAG image:

Back-project the reconstructed 3D catheter path to the left and right CAG imaging planes to obtain the corresponding 2D path. For each point on the 2D path, find the two maxima of the gray gradient in the direction perpendicular to the path value, get the left and right edges of the vessel lumen, and then complete the three-dimensional reconstruction of the entire vessel lumen under the assumption that the cross-section of the lumen is an ellipse, and the result is used to determine the spatial direction of each frame of ultrasound images;

c. Extraction of vessel wall contours in intravascular ultrasound image sequences:

Manually select several feature points on the contour of the inner and outer membranes of the vessel wall in the first frame image, and use the polygon formed by connecting these points as the initial position, obtain the outline of the inner and outer membranes of the vessel wall through snake deformation, and segment the vessel wall and possible plaques, for subsequent frames, the extraction result of the previous frame is used as the initial position of the snake to complete the segmentation of continuous multi-frame IVUS images;

d. Determine the axial position of each frame of the intravascular ultrasound image:

According to the collection order and spacing of IVUS images, arrange the IVUS images of each frame in order along the reconstructed 3D catheter retraction path, and determine the axial position of each frame of images;

e. Determine the spatial orientation of each frame of the intravascular ultrasound image:

The local coordinate system of each frame of ultrasound images is established on the reconstructed 3D catheter path, that is, the Frenet-Serret frame. The three coordinate axes are the unit tangent vector t, the unit principal normal vector n and the unit secondary normal vector b. The position is at the center of the IVUS image;

Rotate each frame of the ultrasound image around the catheter to its correct orientation to determine the spatial orientation of the IVUS image:

Use ρ to represent the centrifugal vector where the center of gravity of the vessel wall contour deviates from the catheter position, project the elliptical contour of the vessel lumen reconstructed from the CAG image onto the corresponding ultrasound image, and use μ to represent the deviation from the centerline of the vessel lumen ellipse contour from the catheter position The centrifugal vector of , ε is the modulus of vector ρ, θ is the angle between ρ and μ, the spatial orientation of the ultrasonic image sequence is determined by statistical optimization method, the purpose is to minimize θ:

Set a moving window with a fixed width w, and perform statistical analysis in this window. For an ultrasound image sequence composed of N frames, there are n _w =N-(w-1) moving windows, and at each window position m, The cumulative eccentric distance ∑ε _m , the average weighted eccentricity angle θ _m and the weighted standard deviation σ(θ _m ) of the eccentricity angle can be calculated by the following formulas:

$\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}}\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&sigma;</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}}\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

Using these values, at each window position, calculate the reliability weighting factor: r _m =∑ε _m /σ(ε _m ), give a larger weighting factor at the position with a larger eccentricity distance, and limit σ(θ _m ) larger position, a corrected eccentricity angle θ _corr is calculated by the following formula:

${\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; corr</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;r</span>}}\underset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; w</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&Sigma;r</span>}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; w</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}$

And apply it to all images of the sequence to obtain the spatial orientation of each frame image;

f. Using the surface extraction method based on NURBS surface fitting to complete the drawing of the inner and outer surfaces of the blood vessel lumen.

The implementation method of the above-mentioned interactive coronary virtual angioscope, the specific steps of using the virtual reality modeling language to realize the visualization of the endoscopic roaming mode of the coronary artery reconstruction results are as follows:

①Calculation of roaming path:

According to the three-dimensional reconstruction results of the ultrasound catheter retraction path, the coordinate P _i of the IVUS image acquisition point of the i-th frame along the catheter is the current viewpoint position, P _i+1 is the next viewpoint position, and the direction vector of the negative semi-axis of the z-axis $-\stackrel{\&Right\; Arrow;}{z}e=\left[\begin{array}{ccc}0& 0& -1\end{array}\right]$ is the initial viewpoint direction, then in VRML, the position vector of the current viewpoint is ${\stackrel{\&Right\; Arrow;}{p}}_i=\left[\begin{array}{ccc}{x}_i& {\mathrm{the\; y}}_i& z_i\end{array}\right],$ The rotation axis is the current viewpoint direction is ${\stackrel{\&Right\; Arrow;}{d}}_i={\stackrel{\&Right\; Arrow;}{p}}_{i+1}-{\stackrel{\&Right\; Arrow;}{p}}_i=\left[\begin{array}{ccc}{d}_{\mathrm{xi}}& d_{\mathrm{yi}}& d_{\mathrm{the\; zi}}\end{array}\right]$ The unit vector of $-{\stackrel{\&Right\; Arrow;}{z}}_i={\stackrel{\&Right\; Arrow;}{d}}_i/\left|\right|{\stackrel{\&Right\; Arrow;}{d}}_i\left|\right|,$ The axis of rotation is ${\stackrel{\&Right\; Arrow;}{r}}_i=\frac{-{\stackrel{\&Right\; Arrow;}{z}}_e\&times;-{\stackrel{\&Right\; Arrow;}{z}}_i}{\left|\right|-{\stackrel{\&Right\; Arrow;}{z}}_e\&times;-{\stackrel{\&Right\; Arrow;}{z}}_i\left|\right|},$ The rotation angle is

②Insert the IVUS pixel data into the virtual scene and display it in a semi-transparent display mode;

③Develop an interactive graphical user interface.

The realization method of the above-mentioned interactive coronary virtual angioscope, the three-dimensional reconstruction method of the ultrasonic catheter retraction path is: first establish the perspective projection imaging model of the CAG system at two approximately vertical angles, and then according to the synchronously recorded during the angiography process The distance and angle parameters are used to obtain the geometric transformation matrix of the imaging system, and then the three-dimensional snake model technology is used to directly deform the snake in space to complete the three-dimensional reconstruction of the catheter retraction path;

In the implementation method of the above-mentioned interactive coronary virtual angioscope, the value range of the included angle between the two acquisition angles of the coronary angiography image is 60° to 120°, and only one shot is taken at the starting point of the ultrasonic catheter retraction path. For contrast images.

The present invention combines the geometrical information of blood vessel space obtained from two single-plane angiography images with approximately orthogonal angles and the lumen cross-section information obtained from intravascular ultrasound images, fully utilizes the complementarity of the two imaging methods, and completes the accurate detection of blood vessels. Three-dimensional reconstruction and visualization of coronary arteries in endoscopic roaming mode using virtual reality modeling language. The invention realizes the interactive visit and display of the three-dimensional blood vessel model, and provides an ideal solution for the development of coronary atherosclerotic lesions, the visual diagnosis and treatment of coronary heart disease, the evaluation of interventional treatment effects, and the training of medical personnel. platform.

Description of drawings

The present invention will be described in further detail below in conjunction with the accompanying drawings.

Fig. 1 is the flow chart of three-dimensional reconstruction blood vessel according to the method of the present invention;

Fig. 2 is a schematic diagram of CAG and IVUS image acquisition according to the method of the present invention;

Fig. 3 is a schematic diagram of imaging at two angles of the imaging system according to the method of the present invention;

Fig. 4 is a schematic diagram of determining the relative orientation of each frame of ultrasonic images according to the method of the present invention;

Fig. 5 is a schematic diagram of the eccentric distance and eccentric angle of the ultrasonic image according to the method of the present invention;

Fig. 6 is a schematic diagram of determining the roaming viewpoint position according to the method of the present invention;

Fig. 7 is a schematic diagram of determining the direction of the roaming viewpoint according to the method of the present invention.

旋转轴；φ_i、旋转角；

z轴负半轴的方向矢量；

待求的视点方向；

y轴的单位矢量，也即VRML中默认的向上的方向；ε_i、在x-y平面内的旋转角。The symbols in the figure are: Image A, Image B, imaging plane; s ₁ , s ₂ , the position of the focal point of the X-ray source during the two imaging processes; s ₁ x ₁ y ₁ z ₁ , the space coordinates with s ₁ as the origin system; s ₂ x ₂ y ₂ z ₂ , space coordinate system with s ₂ as the origin; U ₁ V ₁ O ₁ , rectangular coordinate system on imaging plane A; U ₂ V ₂ O ₂ , rectangular coordinate system on imaging plane B Coordinate system; vertical distance from D ₁ , s ₁ to imaging plane A; vertical distance from D ₂ , s ₂ to imaging plane B; P, point on the space vessel; p ₁ , projection of point P on imaging plane A; p ₂ , the projection of point P on the imaging plane B; the abscissa of u ₁ and p ₁ in the coordinate system U ₁ V ₁ O ₁ ; the ordinate of v ₁ and p ₁ in the coordinate system U ₁ V ₁ O ₁ ; the abscissa of u ₂ and p ₂ in the coordinate system U ₂ V ₂ O ₂ ; the ordinate of v ₂ and p ₂ in the coordinate system U ₂ V ₂ O ₂ ; c(s), representing the space of the 3D catheter path Parameter curve; C, the position of the catheter in the ultrasound image, which is also the center of the ultrasound image; O _C , the center of the ellipse cross-sectional contour (that is, when the cross-section of the blood vessel is assumed to be an ellipse, the corresponding vessel centerline in the three-dimensional reconstruction based on CAG position); O _I , the center of gravity of the lumen cross-sectional profile extracted from the ultrasound image; ρ, ρ=O _I -C is the centrifugal vector that the center of gravity of the ultrasound profile deviates from the catheter; μ, μ=O _C -C is the center of the blood vessel The centrifugal vector of the line deviating from the catheter; the angle between θ, ρ and μ; P _i , P _i+1 , viewpoint; position vector;

Rotation axis; φ _i , rotation angle;

The direction vector of the negative semi-axis of the z-axis;

The direction of the viewpoint to be requested;

The unit vector of the y-axis, which is the default upward direction in VRML; ε _i , Rotation angle in the xy plane.

Symbols used in this paper: t, unit tangent vector; n, unit principal normal vector; b, unit secondary normal vector; ε, modulus of vector ρ; w _, moving window width _; ), the weighted standard deviation of the included eccentricity angle; r _m , the reliability weight factor; θ _corr , the corrected eccentricity angle.

Detailed ways

The steps of the present invention are described in detail below in conjunction with accompanying drawings:

(1) Image acquisition:

Acquisition equipment includes a C-arm single-sided X-ray angiography machine and an intravascular ultrasound imager.

Referring to Figure 2, IVUS and CAG imaging are performed simultaneously. Routine puncture through the right femoral artery or the brachial artery of the upper arm for selective coronary angiography. Under the guidance of the X-ray fluoroscopy image, insert the high-frequency ultrasound probe catheter to the distal end of the blood vessel. After the ultrasound probe was connected to the ultrasound imager to remove artifacts, the catheter was retracted at a constant speed and equidistant by motor control. When the probe catheter rotates 360° at 1800 r/min, real-time vascular section images of 30 frames per second are obtained continuously. The commonly used clinical method of allowing the patient to hold his breath during the withdrawal of the catheter is adopted to reduce the influence of respiratory movement. Ultrasound images are acquired by means of ECG gating, thereby reducing the influence of cardiac motion.

Only at the starting point of the catheter retraction path, a pair of contrast images with approximately vertical angles are taken at the corresponding cardiac phases by means of ECG gating. Due to the mechanical ultrasonic catheter probe, the ultrasonic transducer is located at the head end of a flexible shaft, and the shaft rotates in the outer sheath while the sheath is fixed, so the stability of the withdrawal path can be guaranteed. During the imaging process, the imaging angle and the distance from the focus of the X-ray source to the receiving screen were recorded.

(2) Extraction and three-dimensional reconstruction of catheter path and lumen edge in contrast images:

The present invention first establishes the perspective projection imaging model of the CAG system at two approximately vertical angles (accompanying drawing 3). Afterwards, according to the distance and angle parameters recorded synchronously during the imaging process, the geometric transformation matrix of the imaging system is obtained. Then use the 3D snake model technology to complete the 3D reconstruction of the catheter path.

The snake model, also known as the active contour model, is a deformation model technology proposed by Kass et al. in 1987 (Kass M, Witkin A, Terzopoulos T. Snakes: active contour models. International Journal of Computer Vision, 1987 , 1(4):321-331), it has been widely used in the field of image processing in recent years to complete image segmentation, matching and motion tracking.

The specific implementation method is: the initial position of the snake is obtained by manually picking points, that is, a number of sampling points are manually selected on a projection of the catheter (generally, the starting point, the end point and 3 to 4 intermediate points of the retraction path can be selected), and then according to The epipolar constraints get the corresponding points of these points on another projection. Calculate their three-dimensional coordinates from these groups of corresponding points, connect these 3D points with straight line segments, and use the obtained polyline as the initial position of the 3D snake.

The energy function of the snake model is:

$\mathrm{<span\; class="notranslate">\; E.</span>}<span\; class="notranslate">\; =</span>{<span\; class="notranslate">\; \&Integral;</span>}_{\mathrm{<span\; class="notranslate">\; 0</span>}}^{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; [</span>{\mathrm{<span\; class="notranslate">\; E.</span>}}_{\mathrm{<span\; class="notranslate">\; int</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; E.</span>}}_{\mathrm{<span\; class="notranslate">\; ext</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; the\; s</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; ]</span>\mathrm{<span\; class="notranslate">\; ds</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

where c(s)=(x(s), y(s), z(s)), s∈[0,1] is a cubic B-spline curve representing the catheter. The expression of internal energy E _int in formula (1) is:

E _int (c(s))＝(α|c′(s)| ² +β|c”(s)| ² )/2 (2)

Where c'(s) and c"(s) are the first and second derivatives of c(s) respectively. The internal energy ensures the continuity and smoothness of the curve.

The external energy function E _ext is the external force to ensure the convergence of the snake. It consists of two parts, corresponding to the left and right projections respectively, to ensure that the projection of the three-dimensional curve on the two-angle imaging plane is exactly located at the corresponding catheter projection:

${\mathrm{E.}}_{\mathrm{ext}}=\mathrm{\&gamma;}(I_L({\mathrm{<figure-callout\; id="1"\; label="u"\; filenames="H2008100800020E00032.png"\; state="\{\{state\}\}">u</figure-callout>}}_1,v_1)+I_R({\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="H2008100800020E00031.png"\; state="\{\{state\}\}">u</figure-callout>}}_2,v_2))+\mathrm{\&lambda;}(|\&dtri;I_L({\mathrm{<figure-callout\; id="1"\; label="u"\; filenames="H2008100800020E00032.png"\; state="\{\{state\}\}">u</figure-callout>}}_1,v_1)|+|\&dtri;I_R({\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="H2008100800020E00031.png"\; state="\{\{state\}\}">u</figure-callout>}}_2,v_2)|)$ (3)

$<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&gamma;</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; L</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; L</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&lambda;</span>}<span\; class="notranslate">\; (</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&dtri;</span>{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; L</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; L</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; +</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; \&dtri;</span>{\mathrm{<span\; class="notranslate">\; I</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; f</span>}}_{\mathrm{<span\; class="notranslate">\; R</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; )</span>$

IL(u₁，v₁)分别是左投影点的灰度和灰度梯度值；I_R(u₂，v₂)和

I_R(u₂，v₂)分别是右投影点的灰度和灰度梯度。由于造影图像中，血管的灰度值比背景小，所以权重系数γ取正值。根据透视投影成像的几何关系和外极线约束关系，可推导出u₁、v₁、u₂和v₂都是空间点三维坐标c＝(x₁，y₁，z₁)的函数：where I _L (u ₁ , v ₁ ) and

IL(u ₁ , v ₁ ) are the grayscale and grayscale gradient value of the left projected point respectively; I _R (u ₂ , v ₂ ) and

I _R (u ₂ , v ₂ ) are the grayscale and grayscale gradient of the right projected point, respectively. Since the gray value of the blood vessel is smaller than the background in the contrast image, the weight coefficient γ takes a positive value. According to the geometric relationship of perspective projection imaging and the constraint relationship of epipolar lines, it can be deduced that u ₁ , v ₁ , u ₂ and v ₂ are all functions of the three-dimensional coordinates of space points c=(x ₁ , y ₁ , z ₁ ):

[u ₁ v ₁ ] ^T = F _L (c), [u ₂ v ₂ ] ^T = F _R (c) (4)

Then, by minimizing the energy function of equation (1), the final position of the snake curve determines the three-dimensional axis of the catheter. This method avoids the point-by-point matching between two angles based on the outer pole constraint, and improves the reconstruction accuracy and operation speed.

According to the geometric transformation matrix of the imaging system, the reconstructed 3D catheter path is back-projected to the left and right CAG imaging planes to obtain the corresponding 2D path. For each point on the 2D path, the left and right edges of the vessel lumen are extracted by looking for two maximum values of the gray gradient in the direction perpendicular to the path. Afterwards, on the premise that the cross section of the lumen is assumed to be an ellipse, the three-dimensional reconstruction of the entire lumen is completed, and the result is used for subsequent determination of the spatial direction of each frame of ultrasound images.

(3) Extraction of the vessel wall contour in the intravascular ultrasound image sequence:

The present invention uses a snake model combined with dynamic programming to complete the extraction of the contour of the inner and outer membranes of the blood vessel wall in the IVUS image sequence. The operator only needs to manually select several feature points on the target contour in the first frame, and the polygon formed by connecting these points is used as the initial position of the snake. For subsequent frames, the extraction result of the previous frame is used as the initial position of the snake in the next frame to complete the segmentation of continuous multi-frame IVUS images, which can greatly save calculation time.

(4) Fusion of IVUS and CAG data:

There are two main problems to be solved here: determine the 3D axial position and spatial orientation of each IVUS frame.

(4.1) Determination of the three-dimensional axial position of the ultrasound image:

During the process of acquiring ultrasound images, the catheter is continuously pulled out from the distal end to the proximal end at a uniform speed and at equal intervals by means of motor drive. By adjusting the speed at which the catheter is pulled out, the distance between the cut planes can be adjusted as required. After the axis of the catheter is reconstructed from the CAG image, the axial position of each frame of the IVUS image can be determined by arranging each frame of the IVUS image sequentially in the axial direction according to the known section spacing.

(4.2) Determination of the spatial orientation of each frame of ultrasound images

The present invention uses a non-iterative statistical optimization method to calculate the spatial orientation of each frame of ultrasonic images. First, the local coordinate system of each frame of ultrasound images is established on the reconstructed 3D catheter path, that is, the Frenet-Serret frame, and the three coordinate axes are the unit tangent vector t, the unit principal normal vector n and the unit secondary normal vector b (attached Figure 4), the coordinate origin is the position of the catheter in the IVUS image. After completing the three-dimensional reconstruction of the catheter path, its 3D curve equation c(s) can be obtained. According to the knowledge of differential geometry, t, n and b can be calculated according to the curve equation as follows:

$\left\{\begin{array}{c}\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}\\ \mathrm{<span\; class="notranslate">\; b</span>}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&times;</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}\\ \mathrm{<span\; class="notranslate">\; no</span>}<span\; class="notranslate">\; =</span>\frac{{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span>}{<span\; class="notranslate">\; |</span>{<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; c</span>}<span\; class="notranslate">\; \&prime;</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; the\; s</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; )</span><span\; class="notranslate">\; |</span>}\end{array}\right.<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 5</span>}<span\; class="notranslate">\; )</span>$

Among them, “×” represents the cross product of vectors; “·” represents the point product of vectors; c′(s) and c″(s) are the first-order and second-order derivatives of c(s) respectively.

The position of the catheter is located in the center of the IVUS image, and the center of gravity of the segmented target contour generally does not coincide with the position of the catheter, as shown in Figure 5, where point C represents the catheter, and O _C is the center of the ellipse contour (that is, in the hypothetical vessel cross-section When is an ellipse, the position of the centerline of the blood vessel in the three-dimensional reconstruction based on CAG), 0 ₁ is the center of gravity of the blood vessel cross-sectional contour extracted from the ultrasound image. The degree to which the center of gravity of the profile deviates from the position of the catheter is represented by the centrifugal vector ρ: ρ=O _I −C.

Since the centerline of the blood vessel does not coincide with the path of the catheter, the contour of the ultrasound image at the same position of the blood vessel is inconsistent with the orientation of the ellipse contour reconstructed based on the contrast image, and the ellipse contour is projected onto the corresponding ultrasound image. The centrifugal vector μ of the ellipse contour is also used to represent the degree of deviation of the centerline of the blood vessel from the position of the catheter: μ=O _C −C.

The matching error of the ultrasonic image can be expressed by the modulus ε of the vector ρ and the angle θ between ρ and μ. The present invention uses a statistical optimization method to determine the absolute orientation of the ultrasonic image sequence, with the aim of minimizing the angle θ between the elliptical contour and the centrifugal vector of the ultrasonic contour. Set a moving window with a fixed width w, in which the statistical analysis is performed. For an ultrasound image sequence consisting of N frames, there are n _w =N-(w-1) moving windows. At each window position m, the cumulative eccentric distance ∑ε _m , the average weighted eccentricity angle θ _m and the weighted standard deviation σ(θ _m ) of the eccentricity angle can be calculated by the following formulas:

$\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}}\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&sigma;</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}}\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

Using these values, at each window position, a reliability weighting factor is calculated: r _m =Σε _m /σ(ε _m ). A larger weight factor is given to a position with a larger eccentricity distance, and a position with a larger σ(θ _m ) is restricted at the same time. A corrected eccentricity angle θ _corr is calculated by the following formula:

${\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; corr</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;r</span>}}\underset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; w</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; \&Sigma;r</span>}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; w</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}$

And apply it to all images in the sequence to obtain the spatial orientation of each frame image.

(5) Fitting of the inner and outer surfaces of the vascular lumen

After completing the edge extraction on the IVUS image sequence and determining the spatial position of each frame, the present invention adopts NURBS (Non-Uniform Rational B-spline) surface to fit the sampling points on each cross-section correctly arranged along the three-dimensional retraction path to obtain continuous 3D vessel surface.

(6) Interactive coronary virtual endoscopy system

Using virtual reality modeling language to display the reconstruction results of the endoscopic roaming mode can not only display the overall appearance of the reconstructed vessel segment, but also display the long-axis longitudinal section image.

Including the calculation of the roaming path, the drawing of the reconstructed inner and outer surfaces of the lumen, the display of IVUS image data in the virtual scene and the development of an interactive graphical user interface.

(6.1) Calculation of roaming path

A walk-through path is a sequence of viewpoints within the lumen of the target vessel. For each viewpoint, its position and orientation need to be determined, where the orientation is represented by the rotation of the observer's local coordinate system around an arbitrary axis.

来表示，其中

是位置矢量，

是旋转轴，φ_i是旋转角(如附图6所示)。Determination of the position of the viewpoint: In VRML, the position and direction of a point P _i (namely the viewpoint) on the roaming path use a ternary number

to represent, among them

is the position vector,

is the axis of rotation, and φ _i is the angle of rotation (as shown in Figure 6).

就是沿导管的第i帧IVUS图像采集点的坐标，也即该帧图像的中心坐标；对于后者，就是从第i帧IVUS图像中分割出的管腔轮廓的重心坐标。由于导管的刚性和连续性都比计算出的管腔中心线要好，因此本发明采用第一种方法确定

从而得到更为连续光滑的动画效果。The roaming of the virtual observer in the target vessel can be carried out either along the withdrawal path of the IVUS catheter or along the centerline of the lumen. so ${\stackrel{\&Right\; Arrow;}{p}}_i=\left[\begin{array}{ccc}{x}_i& {\mathrm{the\; y}}_i& z_i\end{array}\right]$ It can be directly obtained according to the 3D reconstruction results of the two: for the former,

is the coordinates of the i-th frame IVUS image collection point along the catheter, that is, the center coordinates of the frame image; for the latter, is the coordinates of the center of gravity of the lumen contour segmented from the IVUS image of the i-th frame. Since the rigidity and continuity of the catheter are better than the calculated centerline of the lumen, the present invention uses the first method to determine

This results in a more continuous and smooth animation effect.

与

的叉积即为旋转轴：Determination of the direction of the viewpoint: the initial value of the direction of the viewpoint is set as the negative half axis along the z axis, as shown in Figure 6, where $-\stackrel{\&Right\; Arrow;}{z}e=\left[\begin{array}{ccc}0& 0& -1\end{array}\right]$ Indicates the direction vector of the negative semi-axis of the z-axis, that is, the direction of the initial viewpoint. Since roaming is carried out along the retraction path of the catheter, according to the method for determining the position of the viewpoint above, the position vector of the next viewpoint P _i+1 can be known from the current viewpoint P _i ${\stackrel{\&Right\; Arrow;}{p}}_{i+1}=\left[\begin{array}{ccc}{x}_{i+1}& {\mathrm{the\; y}}_{i+1}& z_{i+1}\end{array}\right],$ so that the vector ${\stackrel{\&Right\; Arrow;}{d}}_i={\stackrel{\&Right\; Arrow;}{p}}_{i+1}-{\stackrel{\&Right\; Arrow;}{p}}_i=\left[\begin{array}{ccc}{d}_{\mathrm{xi}}& d_{\mathrm{yi}}& d_{\mathrm{the\; zi}}\end{array}\right],$ Its unit vector is the direction of the viewpoint to be sought $-{\stackrel{\&Right\; Arrow;}{z}}_i={\stackrel{\&Right\; Arrow;}{d}}_i/\left|\right|{\stackrel{\&Right\; Arrow;}{d}}_i\left|\right|.$

and

The cross product of is the axis of rotation:

${\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; r</span>}}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>\frac{<span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; z</span>}}}_{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; z</span>}}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span><span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; z</span>}}}_{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; \&times;</span><span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; z</span>}}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; |</span><span\; class="notranslate">\; |</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 6</span>}<span\; class="notranslate">\; )</span>$

与所决定的平面，如附图7所示，视点的初始方向

绕旋转轴

旋转φ_i角，即得到当前的视点方向

即由当前视点P_i指向下一个视点P_i+1的方向，旋转角为：Right now perpendicular to

and The determined plane, as shown in Figure 7, the initial direction of the viewpoint

around the axis of rotation

Rotate the φ _i angle to get the current viewpoint direction

That is, from the current viewpoint P _i to the direction of the next viewpoint P _i+1 , the rotation angle is:

的z分量为0，表示在x-y平面内，所以

也可由y轴的单位矢量

在x-y平面内旋转ε_i得到：According to the calculation formula of vector cross product,

The z component of is 0, which means in the xy plane, so

It can also be represented by the unit vector of the y-axis

Rotating ε _i in the xy plane gives:

${\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; r</span>}}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>{\stackrel{<span\; class="notranslate">\; \&Right\; Arrow;</span>}{\mathrm{<span\; class="notranslate">\; the\; y</span>}}}_{\mathrm{<span\; class="notranslate">\; e</span>}}{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; z</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 8</span>}<span\; class="notranslate">\; )</span>$

是VRML中默认的向上的方向。where R _z (ε _i ) represents the rotation matrix around the z axis,

is the default up direction in VRML.

At the end of the roaming path, that is, for a point P _i in the set of viewpoints {P ₀ , P ₁ ,...,P _n-1 }, when i=n-1, since there is no P _i+1 , the aforementioned method to calculate the rotation axis and rotation angle, at this time: ${\stackrel{\&Right\; Arrow;}{r}}_{\mathrm{no}-1}={\stackrel{\&Right\; Arrow;}{r}}_{\mathrm{no}-2},$

(6.2) Display of IVUS pixel data in VRML

The present invention adopts a semi-transparent display mode for the IVUS pixel data inserted in the virtual scene, that is, the transparency value of each pixel in a frame of IVUS image is not the same constant, but depends on the position of the pixel in the image and its gray value. In the ultrasound image, except for the vessel wall and plaque, other structures should be invisible in the virtual endoscopic scene. At this time, using the previous two-dimensional segmentation results of the ultrasound image, the echo signal outside the lumen and adventitia will be represented The pixels of are set to be fully transparent, allowing the walk path to traverse these areas unobstructed. For regions of interest such as tube walls and plaques, the transparency value depends on the gray value of the pixel, for example, a bright echo signal indicates a possible plaque, so its opacity value should be set to a higher value; in the image Dark areas in may represent lumens of other vessels or other structures that are not echogenic, and should be set to a lower value for their opacity values.

(6.3) Interactive graphical user interface

The present invention designs and develops a concise, clear, convenient and flexible user control panel in the VRML environment, so that it can complete the following functions: ① The user can open and close the control panel at any time, and when opening, the target scene should not be blocked as much as possible. ②When the virtual observer advances along the roaming path, at a certain point of view, the user can switch between different display modes, for example: display the IVUS image acquired at this point in the correct direction and position; or only display this point The luminal surface at the location (the semi-transparent IVUS frame can be turned on or off at the same time); or the luminal surface that has been color-coded, where the color coding indicates quantitative measurement results, etc. At the same time, the user can enter or exit the observation mode of the virtual endoscope at any time to display the overall appearance of the reconstructed vessel segment, or the long-axis longitudinal section image. ③ The roaming speed and direction can be adjusted, and the virtual observer can stay at any position in the lumen.

Claims (3)

1. A method for realizing an interactive coronary virtual angioscope, characterized in that it combines the lumen three-dimensional geometry information obtained by approximately orthogonal X-ray coronary angiography CAG images with the lumen obtained by intravascular ultrasound The cross-sectional data are fused to obtain a three-dimensional model of the blood vessel, and then the virtual reality modeling language VRML is used to interactively describe the blood vessel model to realize the visualization of the coronary artery in the endoscopic roaming mode. The specific steps are as follows: A. Using simultaneously acquired intravascular ultrasound IVUS images and X-ray coronary angiography CAG images to establish a three-dimensional model of blood vessels; The intravascular ultrasound IVUS image is obtained by placing the mechanical ultrasound catheter probe at the distal end of the vessel segment of interest and withdrawing the guide wire at a constant speed and equidistantly. Equidistant intravascular ultrasound IVUS image sequence acquired at the same heart phase, the X-ray coronary angiography CAG image is recorded at the starting point of the catheter retraction path by using a C-arm single-sided X-ray angiography machine Two X-ray coronary angiography CAG images at approximately vertical angles; The specific steps for establishing a three-dimensional model of blood vessels are as follows: a. Based on two X-ray coronary angiography CAG images at approximately perpendicular angles, a three-dimensional ultrasound catheter withdrawal path is reconstructed; b. Reconstruction of three-dimensional vessel lumen from X-ray coronary angiography CAG images: Back-project the reconstructed three-dimensional ultrasound catheter withdrawal path to the imaging planes of the left and right X-ray coronary angiography CAG images to obtain the corresponding two-dimensional ultrasound catheter withdrawal path. For each point, by looking for the two maximum values of the gray gradient in the direction perpendicular to the retraction path of the two-dimensional ultrasound catheter, the left and right edges of the vascular lumen are obtained, and then on the premise that the cross section of the lumen is an ellipse, complete Reconstruction of the entire three-dimensional vessel lumen, the three-dimensional vessel lumen is used to subsequently determine the spatial direction of each frame of intravascular ultrasound IVUS images; c. Extraction of vessel wall contours in IVUS image sequences: In the first frame of intravascular ultrasound IVUS image, manually select several feature points on the contour of the inner and outer membranes of the vessel wall, and use the polygon formed by connecting these feature points as the initial position, and obtain the outline of the inner and outer membranes of the vessel wall through snake deformation , to segment the vessel wall and possible plaques, and for subsequent frames, the extraction result of the previous frame is used as the initial position of the snake to complete the segmentation of continuous multi-frame intravascular ultrasound IVUS images; d. Determine the axial position of each frame of intravascular ultrasound IVUS image: According to the acquisition sequence and spacing of intravascular ultrasound IVUS images, arrange each frame of intravascular ultrasound IVUS images sequentially along the reconstructed three-dimensional ultrasound catheter withdrawal path, and determine the axial position of each frame of intravascular ultrasound IVUS images; e. Determine the spatial orientation of each frame of intravascular ultrasound IVUS image: The local coordinate system of each frame of intravascular ultrasound IVUS images is established on the reconstructed three-dimensional ultrasound catheter retraction path, that is, the Frenet-Serret frame, and the three coordinate axes are the unit tangent vector t, the unit principal normal vector n and the unit vice Normal vector b, the position of the catheter is located in the center of the IVUS image; Rotate each frame of the IVUS image around the catheter to its correct orientation to determine the spatial orientation of the IVUS image: Use ρ to represent the centrifugal vector where the center of gravity of the vessel wall contour deviates from the position of the catheter, project the elliptical contour of the vessel lumen reconstructed from the X-ray coronary angiography CAG image onto the corresponding intravascular ultrasound IVUS image, and use μ to represent the vessel The center line of the lumen ellipse outline deviates from the centrifugal vector of the catheter position, ε is the modulus of the vector ρ, and θ is the angle between ρ and μ. The spatial orientation of the intravascular ultrasound IVUS image sequence is determined by the statistical optimization method, and the purpose is to minimize θ: Set a moving window with a fixed width w, and perform statistical analysis in this window. For the intravascular ultrasound IVUS image sequence composed of N frames, there are n _w =N-(w-1) moving windows, and at each window position At m, the cumulative eccentric distance ∑ε _m and the average weighted eccentric angle and the weighted standard deviation σ(θ _m ) of the eccentric angle are calculated by the following formula: $\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>$ ${\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}}\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>$ $\mathrm{<span\; class="notranslate">\; \&sigma;</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}}\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; m</span>}}{\overset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; w</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; \&epsiv;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$ Using these values, at each window position, calculate the reliability weighting factor: r _m =∑ε _m /σ(ε _m ), give a larger weighting factor at the position with a larger eccentricity distance, and limit σ(θ _m ) larger position, a corrected eccentricity angle is calculated by the following formula ${\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; corr</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;r</span>}}\underset{\mathrm{<span\; class="notranslate">\; m</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; w</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; m</span>}}{\stackrel{<span\; class="notranslate">\; \&OverBar;</span>}{\mathrm{<span\; class="notranslate">\; \&theta;</span>}}}_{\mathrm{<span\; class="notranslate">\; m</span>}}<span\; class="notranslate">\; ,</span>$ $\mathrm{<span\; class="notranslate">\; \&Sigma;r</span>}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 0</span>}}{\overset{{\mathrm{<span\; class="notranslate">\; no</span>}}_{\mathrm{<span\; class="notranslate">\; w</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}}}{\mathrm{<span\; class="notranslate">\; r</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}$ And apply it to all the intravascular ultrasound IVUS images of the sequence, so as to obtain the spatial orientation of each frame of intravascular ultrasound IVUS images; f. Using the surface extraction method based on NURBS surface fitting to complete the drawing of the blood vessel surface; B. Using the virtual reality modeling language VRML to realize the visualization of the endoscopic roaming mode of the coronary artery reconstruction results: Specific steps are as follows: ①Calculation of roaming path: According to the reconstruction result of the three-dimensional ultrasound catheter retraction path, the coordinate Pi of the intravascular ultrasound IVUS image acquisition point along the i-th frame of the catheter is the current viewpoint position, P _i+1 is the next viewpoint position, and the negative semi-axis of the z-axis is direction vector is the initial viewpoint direction, then in VRML, the position vector of the current viewpoint is

The rotation axis is the current viewpoint direction is

The unit vector of The axis of rotation is

The rotation angle is

②Insert the IVUS image pixel data into the virtual scene and display it in a translucent display mode; ③Develop an interactive graphical user interface. 2. according to the realization method of the described interactive coronary virtual angioscope of claim 1, it is characterized in that, the reconstruction method of described three-dimensional ultrasonic catheter retraction path is: at first establish X-ray coronary angiography CAG image system in two approximate Perspective projection imaging model of vertical angle, and then according to the distance and angle parameters recorded synchronously during the imaging process, the geometric transformation matrix of the imaging system is obtained, and then the three-dimensional snake model technology is used to directly deform the snake in space to complete the retraction of the three-dimensional ultrasound catheter Path reconstruction. 3. The method for realizing the interactive coronary virtual angioscope according to claim 2, wherein the value range of the included angle between the two acquisition angles of the X-ray coronary angiography CAG image is 60° to 120° °, and only a pair of contrast images were taken at the starting point of the three-dimensional ultrasound catheter retraction path.

Images