ICG fluorescence image and MRCP three-dimensional model registration fusion display method and system

Abstract

The invention relates to the technical field of medical image information processing, discloses an ICG fluorescence image and MRCP three-dimensional model registration fusion display method and system, effectively solves the problems of deep biliary tract structure blurring and form distortion caused by liver tissue scattering effect in intraoperative fluorescence imaging, uses a preoperative model as anatomical prior, and improves the accuracy of the intraoperative fluorescence imaging. A viewpoint correlation optical path depth field is constructed through reverse ray tracing, a spatial variation point spread function is generated accordingly, physical reverse recovery of scattering halo in a fluorescence image is achieved, non-specific background noise such as vascular leakage is effectively inhibited in combination with topological gravitational constraint of a central skeleton, and the detection accuracy is improved. The real pipe diameter and the center line position of the deep bile duct are accurately recovered, non-rigid registration correction is performed on the three-dimensional model based on the recovered high-fidelity image, and projection errors caused by respiratory movement and soft tissue deformation are eliminated.

Description

Technical Field

[0001] This application relates to the field of medical image information processing technology, and in particular to a method and system for registering and fusion displaying ICG fluorescence images and MRCP three-dimensional models. Background Technology

[0002] Hepatobiliary surgery, due to its complex anatomy and frequent vascular and bile duct variations, has long been considered a high-risk surgical procedure. With the deep integration of digital medicine and minimally invasive techniques, 3D reconstruction technology based on medical imaging data has been widely used in preoperative planning. Among these, magnetic resonance cholangiopancreatography (MRI) provides high-contrast 3D structures of the biliary tree non-invasively, becoming the gold standard for assessing anatomical variations. During the intraoperative phase, indocyanine green fluorescence imaging, with its excellent tissue penetration and specific metabolic characteristics in the near-infrared band, can visualize biliary structures in real time, becoming an important tool to assist surgeons in identifying anatomical landmarks and preventing iatrogenic biliary tract injuries. To further overcome the limitations of intraoperative visualization, current advanced navigation solutions tend to register and fuse high-precision preoperative 3D MRI models with real-time intraoperative 2D fluorescence images. Augmented reality technology projects deep bile ducts or minute variations that cannot be directly observed onto the laparoscopic monitor, providing surgeons with intuitive anatomical guidance.

[0003] However, in actual clinical applications, the liver parenchyma, as a highly scattering and turbid medium rich in mitochondria and heme, exhibits a strong scattering effect on near-infrared photons. When the target bile duct is located at a certain depth below the liver capsule, the fluorescent photons accumulated within the duct undergo complex multiple random scatterings as they penetrate the upper liver tissue to reach the imaging sensor. This results in the originally clear and sharp tubular signal appearing as a diffuse, blurred halo. This optical phenomenon causes thinner bile ducts to appear as bright spots with expanded boundaries and significantly increased diameters in two-dimensional images. Existing image processing techniques for extracting bile duct features often employ gradient detection based on pixel brightness or adaptive threshold segmentation algorithms. These traditional algorithms struggle to physically distinguish between the true anatomical edge and the halo-like virtual image produced by scattering effects, often mistakenly identifying the diffused halo boundary as the bile duct wall. More seriously, due to variations in liver surface curvature and uneven tissue thickness distribution, this morphological distortion caused by scattering is anisotropic, resulting in a nonlinear spatial offset between the extracted bile duct centerline and the true anatomical axis. This distortion of basic data can directly mislead subsequent non-rigid registration calculations, causing the accurate 3D model to be incorrectly stretched or distorted in order to match the distorted fluorescence features. Ultimately, this results in serious positioning errors in augmented reality fusion display, greatly increasing the risk of misjudging anatomical depth or even damaging critical ducts during surgery. Summary of the Invention

[0004] This application proposes a method and system for registering and fusing ICG fluorescence images with MRCP three-dimensional models to solve the problems mentioned in the background art.

[0005] To achieve the above objectives, this application adopts the following technical solution: a method for registering and fusing ICG fluorescence images with MRCP three-dimensional models, comprising the following steps:

[0006] Step S1: Simultaneously acquire the preoperative reconstructed MRCP three-dimensional biliary tract model, the intraoperative real-time acquired two-dimensional ICG fluorescence image and the corresponding laparoscopic camera pose matrix. Based on the pose matrix, perform viewpoint projection on the MRCP three-dimensional biliary tract model. Calculate the effective scattering optical path of the anatomical structure corresponding to each pixel on the imaging plane in the liver tissue through reverse ray tracing logic, and generate a viewpoint-related optical path depth mapping map containing tissue thickness information.

[0007] Step S2: Using the viewpoint-related optical path depth mapping map generated in step S1 as an index, query the preset depth-scattering response model, independently map the corresponding point spread function for each pixel on the imaging plane, and construct a spatially variable point spread function array in which the diffusion scale parameter changes nonlinearly with the effective scattering optical path.

[0008] Step S3: Call the spatial variation point diffusion function array constructed in step S2 as the inverse restoration operator, and use the central skeleton projection line of the MRCP three-dimensional biliary tract model in step S1 as the topological gravity constraint to perform restricted inverse deconvolution processing on the two-dimensional ICG fluorescence image obtained in step S1, eliminate the scattered halo and generate the restored biliary tract image.

[0009] Step S4: Extract the anatomical contour features of the reconstructed bile duct image generated in step S3, correct the spatial registration relationship between the MRCP three-dimensional bile duct model and the reconstructed bile duct image based on the anatomical contour features, and render and overlay the MRCP three-dimensional bile duct model onto the reconstructed bile duct image.

[0010] Furthermore, in step S1, the specific execution process of the reverse ray tracing logic includes:

[0011] Using the optical center of the laparoscopic camera as the origin, virtual rays are emitted to each pixel on the imaging plane to perform multi-level collision detection on the MRCP three-dimensional biliary model.

[0012] Calculate the first intersection point between the virtual ray and the liver capsule mesh in the MRCP three-dimensional biliary tract model and mark it as the media incident point;

[0013] Calculate the second intersection point of the same ray with the intrahepatic bile duct grid and mark it as the scattering termination point;

[0014] Offset sampling is performed along the path connecting the incident point and the scattering termination point of the medium. Based on the pre-annotated anatomical semantic labels in the MRCP three-dimensional biliary model, non-scattering medium regions belonging to portal vein vessels or liver cysts are identified along the path. The length of the non-scattering medium region is subtracted from the Euclidean distance between the incident point and the scattering termination point to obtain the pure geometric penetration depth containing only highly scattered liver parenchyma components.

[0015] Furthermore, in step S1, during the process of generating the viewpoint-related optical path depth map, the calculation of the effective scattered optical path includes weighted correction logic based on the patient's physiological heterogeneity:

[0016] Obtain the patient's preoperative magnetic resonance imaging data or body mass index, extract the patient-specific fat scattering coefficient that characterizes liver fat content, and the magnetic resonance signal intensity normalization index that characterizes the density of local liver tissue.

[0017] Physiological scattering enhancement factors were constructed using patient-specific fat scattering coefficients and normalized magnetic resonance signal intensity indices. These factors were then multiplicatively weighted to the obtained pure geometric penetration depth to reflect the different effects of different fat densities on photon scattering behavior.

[0018] The angle between the virtual ray and the normal vector of the liver capsule mesh at the incident point of the medium is calculated. A Fresnel incident efficiency constraint function is introduced to numerically compensate for the light energy loss caused by large-angle grazing. Finally, the value after geometric semantic elimination, physiological density weighting and incident angle correction is determined as the effective scattered optical path.

[0019] Furthermore, in step S2, the depth-scattering response model employs a logistic scattering saturation calculation logic based on physical optics to establish a nonlinear mapping relationship between the effective scattering optical path and the dispersion scale parameter:

[0020] Set the system diffraction limit base value to characterize the inherent optical properties of the endoscope lens, and the maximum scattering saturation threshold to characterize the photon energy depletion of deep tissue;

[0021] The numerical difference between the effective scattered optical path output in step S1 and the preset ballistic light to scattered light turning optical path is calculated, and this numerical difference is used as an input variable and substituted into the exponential decay function. The nonlinear growth ratio is calculated in combination with the preset scattering phase transition rate factor.

[0022] The initial dispersion scale parameter is obtained by multiplying the maximum scattering saturation threshold by the nonlinear growth ratio and then superimposing the product onto the system diffraction-limited basis value.

[0023] The geometric interaction data between the virtual ray output from the reverse ray tracing logic in step S1 and the liver capsule is called to calculate the cosine value of the incident angle and construct anisotropic morphology correction coefficients. The initial dispersion scale parameters are multiplied to correct the elliptical spot effect caused by oblique incidence. The final dispersion scale parameters are output to simulate the phase transition process of photons from the quasi-ballistic state to the multiple scattering state.

[0024] Furthermore, in step S2, the process of constructing the spatial variation point diffusion function array performs an energy decoupling operation between the ballistic light component and the diffuse light component:

[0025] The ballistic beam weights that retain high-frequency edge information and the diffuse beam weights that cause background blur are calculated based on the effective scattered optical path.

[0026] The ballistic light weight is set to decrease exponentially with the increase of the effective scattered optical path, and the remaining energy after subtracting the ballistic light weight from the total energy is allocated to the diffuse light weight.

[0027] When synthesizing the point spread function kernel of a single pixel, the ballistic light core representing the sharp signal is constructed using the unit impulse function, and the long-tailed distribution smooth function generated by the calculated dispersion scale parameter is used to construct the diffuse light periphery representing the scattered halo.

[0028] Then, the two are multiplied by their respective weights and then linearly superimposed.

[0029] Traverse all pixel positions on the imaging plane, and spatially tile and tensor reassemble each generated point spread function kernel with independent size and shape according to pixel coordinates;

[0030] The final result is a spatial variation point spread function array with a four-dimensional structure. Each element in the spatial variation point spread function array corresponds precisely to the optical degradation characteristics of that location in the image after being affected by physical scattering.

[0031] Furthermore, in step S3, the topological gravitational constraint is established using the technique of constructing a dissected potential energy field:

[0032] Extract the central skeleton projection line generated by projecting the MRCP three-dimensional bile duct model through reverse ray tracing logic in step S1, and use it as the standard spatial reference path for anatomy.

[0033] For each pixel in the imaging plane resolution grid, the Euclidean distance from that pixel to the nearest central skeleton projection line segment is calculated, and this Euclidean distance is defined as the skeleton Euclidean deviation; the skeleton Euclidean deviation characterizes the geometric distance of the current fluorescence signal point from the ideal anatomical path;

[0034] A potential energy attenuation weight field is constructed based on the preset anatomical tolerance radius. This anatomical tolerance radius defines the maximum allowable spatial registration error limit between the MRCP three-dimensional biliary tract model and the real biliary tract in the surgical scenario.

[0035] The ratio of the Euclidean deviation of the skeleton to the anatomical tolerance radius is used as input and substituted into a high-order reciprocal polynomial function with the gravitational field decay order as the power to calculate the potential energy decay weight value between zero and one.

[0036] Furthermore, in step S3, the constrained inverse deconvolution process employs a dual-constraint iterative update mechanism based on tensor product:

[0037] In each iteration cycle, the spatial variation point spread function array constructed in step S2 is first called to perform forward physical scattering simulation on the currently estimated restored image. That is, for each pixel position, its own point spread function kernel is used for local integration to generate a simulated scattering image.

[0038] The brightness ratio of the two-dimensional ICG fluorescence image and the simulated scattering image obtained in step S1 is calculated, and the brightness ratio is back-projected back to the source space through the adjoint inverse operator to construct the physical likelihood term;

[0039] Introducing gradient direction consistency constraint logic, the brightness gradient direction vector of the restored image in the current iteration step is calculated, and the cosine similarity of the angle between the brightness gradient direction vector and the normal vector of the central skeleton projection line is calculated to construct the gradient skeleton consistency tensor.

[0040] The physical likelihood term, the potential energy decay weight field, and the gradient skeleton consistency tensor are subjected to joint multiplication. A topological penalty factor is introduced to adjust the intensity of the intervention of the anatomical prior on the physical reconstruction. Through multiple rounds of iterative updates, the fluorescence signal energy is forced to move towards the skeleton in spatial position and conform to the tubular structure characteristics in morphology until the relative entropy change rate between two adjacent iterations is lower than a preset threshold.

[0041] Furthermore, in step S4, the specific logic for performing geometric consistency verification includes the construction process of the local confidence graph:

[0042] A multi-scale Hessian matrix filter is applied to the reconstructed bile duct image generated in step S3 to extract response feature maps that can enhance the contrast of tubular structures and suppress sheet noise. The reconstructed skeleton binary map reflecting the center line of the tubular structure is generated through binarization.

[0043] Based on the real-time pose matrix obtained in step S1, a binary projection skeleton map of the MRCP 3D biliary tract model is generated at the current viewpoint.

[0044] Establish a local sliding window centered on the current pixel, and calculate the normalized cross-correlation coefficient between the restored skeleton binary map and the projected skeleton binary map within the local sliding window to quantify the local similarity of the two in terms of topological morphology.

[0045] The gradient magnitude of the reconstructed bile duct image is calculated and a gradient significance factor is constructed. Multiplication is performed, and the gradient significance factor is jointly weighted with the normalized cross-correlation coefficient to generate a local confidence map. Each value in the local confidence map accurately represents the degree of topological fit between the image features and model features at the corresponding pixel location. Furthermore, by introducing the gradient significance factor, the confidence of flat regions with excessively low image gradient magnitudes is forced to be reduced to zero.

[0046] Furthermore, in step S4, the logic for correcting the spatial registration relationship and rendering overlay adopts a view plane non-rigid deformation compensation and visual blocking rendering mechanism:

[0047] In the local confidence map, high confidence regions with values ​​higher than the preset safety threshold are selected, and the iterative nearest neighbor algorithm is used in the high confidence region to search for the nearest neighbor correspondence between the skeleton point set of the restored biliary tract image and the projection skeleton point set of the MRCP three-dimensional biliary tract model. The planar displacement difference between each pair of feature points is calculated to construct the residual deviation vector field.

[0048] The residual deviation vector field is fitted globally using a thin plate spline transformation model to construct a two-dimensional mapping function that can describe the nonlinear deformation of the entire field. This two-dimensional mapping function is then used to update the screen output coordinates of the mesh vertices of the MRCP three-dimensional biliary tract model in the rendering pipeline.

[0049] In the final rendering stage, a hierarchical display logic is introduced. For areas with a confidence level higher than the safety threshold, the model outline after deformation compensation is depicted with a highlighted solid line. For areas with a confidence level lower than the safety threshold, a pixel discard operation is performed to prevent the rendering of the model outline, and a visual blocking mechanism is triggered to generate a dynamically flashing semi-transparent warning cloud map to cover the area, so as to intuitively remind the operator that there is registration uncertainty here.

[0050] The ICG fluorescence image and MRCP 3D model registration and fusion display system includes: a viewpoint-related optical path depth map generation module, a spatial variation point diffusion function array construction module, a restricted inverse deconvolution restoration module, and a registration and fusion display module, wherein;

[0051] The viewpoint-related optical path depth mapping generation module is configured to simultaneously acquire the preoperatively reconstructed MRCP three-dimensional biliary tract model, the intraoperatively acquired two-dimensional ICG fluorescence image and the corresponding laparoscopic camera pose matrix, perform viewpoint projection on the MRCP three-dimensional biliary tract model based on the pose matrix, and calculate the effective scattering optical path of the anatomical structure corresponding to each pixel on the imaging plane in the liver tissue through reverse ray tracing logic, and generate a viewpoint-related optical path depth mapping containing tissue thickness information.

[0052] The spatial variation point diffusion function array construction module is configured to use the viewpoint-related optical path depth mapping generated by the viewpoint-related optical path depth mapping generation module as an index to query the preset depth-scattering response model, independently map the corresponding point diffusion function for each pixel on the imaging plane, and construct a spatial variation point diffusion function array whose diffusion scale parameter changes nonlinearly with the effective scattering optical path.

[0053] The constrained inverse deconvolution restoration module is configured to call the spatial variation point diffusion function array constructed by the spatial variation point diffusion function array construction module as the inverse restoration operator, and use the central skeleton projection line of the MRCP three-dimensional biliary tract model generated by the viewpoint-related optical path depth map generation module as the topological gravity constraint to perform constrained inverse deconvolution processing on the two-dimensional ICG fluorescence image obtained by the viewpoint-related optical path depth map generation module to eliminate scattered halo and generate restored biliary tract image;

[0054] The registration and fusion display module is configured to extract the anatomical contour features of the restored bile duct image generated by the restricted inverse deconvolution restoration module, correct the spatial registration relationship between the MRCP three-dimensional bile duct model and the restored bile duct image based on the anatomical contour features, and render and overlay the MRCP three-dimensional bile duct model onto the restored bile duct image.

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

[0056] This invention effectively solves the problem of blurred and morphologically distorted deep bile duct structures caused by the scattering effect of liver tissue in intraoperative fluorescence imaging. It utilizes a preoperative model as an anatomical prior, constructs a viewpoint-related optical path depth field through reverse ray tracing, and generates a spatial variation point diffusion function accordingly. This achieves the physical inverse restoration of scattered halos in fluorescence images. Combined with the topological gravity constraint of the central skeleton, it effectively suppresses non-specific background noise such as vascular leakage, accurately restores the true diameter and centerline position of deep bile ducts, and performs non-rigid registration correction on the three-dimensional model based on the restored high-fidelity image. This eliminates projection errors caused by respiratory motion and soft tissue deformation, and introduces a visual blocking mechanism to shield low-confidence areas, significantly improving the geometric fit of the image augmented reality navigation system and the safety of the surgery. Attached Figure Description

[0057] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on the provided drawings without creative effort:

[0058] Figure 1 This is a flowchart of the method of the present invention;

[0059] Figure 2 This is a system framework diagram of the present invention. Detailed Implementation

[0060] 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, and 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.

[0061] Example 1

[0062] like Figure 1 As shown, this invention provides a method for registering and fusing ICG fluorescence images with MRCP three-dimensional models, comprising the following steps:

[0063] Step S1: Simultaneously acquire the preoperatively reconstructed MRCP three-dimensional biliary tract model, the intraoperatively acquired two-dimensional ICG fluorescence image and the corresponding laparoscopic camera pose matrix. Based on the pose matrix, perform viewpoint projection on the MRCP three-dimensional biliary tract model. Calculate the effective scattering optical path of the anatomical structure corresponding to each pixel on the imaging plane in the liver tissue through reverse ray tracing logic, and generate a viewpoint-related optical path depth mapping map containing tissue thickness information.

[0064] Next, the core technical logic of step S1 in the method of this invention will be explained in detail, namely, how to accurately calculate the effective transmission path of photons inside the liver parenchyma from both geometric and physical dimensions.

[0065] Specifically, in step S1, the execution process of the reverse ray tracing logic includes:

[0066] Using the optical center of the laparoscopic camera as the origin, virtual rays are emitted to each pixel on the imaging plane to perform multi-level collision detection on the MRCP three-dimensional biliary tract model.

[0067] Specifically, the system first reads the pose matrix of the laparoscopic camera that is updated in real time, and constructs a virtual view frustum that is consistent with the parameters of the physical camera.

[0068] For each pixel coordinate in the resolution grid of the imaging sensor, the system generates a unit direction vector that originates from the optical center and points to the spatial location of that pixel.

[0069] At the same time, the system calls the preoperatively reconstructed MRCP three-dimensional biliary tract model to ensure that the MRCP three-dimensional biliary tract model has been accurately transformed to the current camera coordinate system through rigid registration, thereby providing a unified spatial reference for subsequent geometric collision calculations.

[0070] First, calculate the first intersection point between the virtual ray and the liver capsule mesh in the MRCP three-dimensional biliary tract model and mark it as the media incident point.

[0071] In this step, the system uses a ray-triangle intersection detection algorithm to traverse the geometric mesh on the surface of the liver capsule.

[0072] The system selects the point with the closest Euclidean distance to the optical center among all intersection points and defines it as the physical refraction interface where photons enter the liver parenchyma from the air medium, i.e., the incident point of the medium.

[0073] The second intersection of the same ray with the intrahepatic bile duct grid is then calculated and marked as the scattering termination point.

[0074] The system controls the virtual ray to continue extending in its original direction after passing through the liver capsule, searching for intersections with the grid on the surface of the intrahepatic bile duct tree inside the liver.

[0075] The system defines the first point of collision between the ray and the bile duct grid as the emission source location of the fluorescence signal, i.e., the scattering termination point; if the ray does not intersect with the bile duct grid, the pixel is marked as an invalid region and no further optical path calculation is performed.

[0076] Based on this, discrete step sampling is performed on the path connecting the incident point and the scattering termination point of the medium. According to the pre-annotated anatomical semantic labels in the MRCP three-dimensional biliary model, the non-scattering medium region belonging to the portal vein or liver cyst on the path is identified, and the length of the non-scattering medium region is subtracted from the Euclidean distance between the incident point and the scattering termination point to obtain the pure geometric penetration depth containing only the highly scattered liver parenchyma components.

[0077] In practice, the system sets an offset step sampling length on the line segment formed by the incident point and the scattering termination point of the medium and performs point-by-point scanning. The offset step sampling length represents the resolution of spatial geometric sampling. In this embodiment, it is set to 0.5 mm. The value is based on the average diameter of the secondary branches of the portal vein to prevent missing small blood vessels and thus ensure that the medium components on the path can be accurately identified.

[0078] The system reads the semantic label of each voxel in the MRCP model. When the sampling point falls within the voxel range of portal vein (manifested as light absorption) or liver cyst (manifested as light transmission), the system includes the step size in the total length of the non-scattering path.

[0079] Finally, the system performs a subtraction operation, subtracting the total accumulated non-scattering path length from the total Euclidean distance to obtain the pure geometric penetration depth, which represents only the thickness of the liver parenchyma.

[0080] This pure geometric penetration depth serves as the basis for subsequent calculations of the effective scattered optical path, ensuring that physical optical path calculations exclude interference from non-uniform media.

[0081] Through the geometric elimination operation described above, the algorithm eliminates the optical path calculation error caused by ignoring the differences in the medium, providing an accurate geometric basis for subsequent optical physics calculations.

[0082] In step S1, during the generation of the viewpoint-related optical path depth map, the calculation of the effective scattered optical path includes weighted correction logic based on patient physiological heterogeneity:

[0083] This is because geometric depth alone cannot reflect the significant differences in the optical scattering properties of liver tissue among different patients, and physiological parameters need to be introduced for dimensional correction.

[0084] First, obtain the patient's preoperative magnetic resonance imaging data or body mass index, extract the patient-specific fat scattering coefficient that characterizes liver fat content, and the normalized index of magnetic resonance signal intensity that characterizes the density of local liver tissue.

[0085] The patient-specific fat scattering coefficient represents the shortening effect of lipid droplets within hepatocytes on the mean free path of near-infrared photon scattering. For patients with severe fatty liver, a value of 0.6 to 0.8 is preferred, while for normal livers, a value of 0.1 to 0.2 is preferred. This value is based on linear regression data of fat content and reduced scattering coefficient from literature on the optical properties of biological tissues. Its purpose is to amplify the optical path value of patients with fatty liver, thus matching a stronger scattering blurring effect.

[0086] The normalized index for magnetic resonance signal intensity represents the relative water content and density of local liver tissue. In this embodiment, its value is normalized to a floating-point number between zero and one. Its function is to act as a local fine-tuning factor, reflecting the subtle influence of tissue edema or fibrosis on scattering.

[0087] Subsequently, a physiological scattering enhancement factor was constructed using the patient-specific fat scattering coefficient and the normalization index of magnetic resonance signal intensity. The pure geometric penetration depth obtained in the above steps was then multiplied and weighted to reflect the different effects of different fat densities on photon scattering behavior.

[0088] The specific operational logic is as follows: The system performs a multiplication operation to calculate the product of the patient-specific fat scattering coefficient and the normalized index of the magnetic resonance signal intensity to obtain the scattering increment value. Then, the scattering increment value is added to the unit reference value to obtain the physiological scattering enhancement factor. Finally, the system multiplies the pure geometric penetration depth with the physiological scattering enhancement factor to obtain the corrected optical path length.

[0089] This step ensures that high-fat-density liver tissue is mapped to a larger effective optical path, thereby triggering a larger deblurring kernel.

[0090] Simultaneously, the angle between the virtual ray and the normal vector of the liver capsule mesh at the incident point of the medium is calculated, and the Fresnel incident efficiency constraint function is introduced to numerically compensate for the light energy loss caused by large-angle grazing. Finally, the value after geometric semantic elimination, physiological density weighting and incident angle correction is determined as the effective scattered optical path.

[0091] In practical implementation, the system calculates the dot product of the ray direction and the normal vector of the coating surface to obtain the cosine value of the incident angle.

[0092] Regarding the incident angle threshold in the Fresnel incident efficiency constraint function, this incident angle threshold is the critical angle that distinguishes effective refracted light from ineffective grazing light (total internal reflection). In this embodiment, it is set to 75 degrees. Its value is based on Fresnel's law of reflection. After exceeding this angle, the light energy entering the tissue will decrease exponentially. Its function is that when the incident angle is greater than this threshold, the function output approaches zero weight, suppressing the false deep signal generated by grazing in the liver edge region.

[0093] The system multiplies the physiologically corrected optical path value with the constraint weight to obtain the final effective scattered optical path, and stores it in the viewpoint-related optical path depth map.

[0094] Finally, through the above series of rigorous logical operations, the system outputs a viewpoint-related optical path depth mapping map that is strictly aligned with the current laparoscopic viewpoint. Each pixel value in this map represents the precise and effective scattered optical path after geometric correction, physiological weighting, and physical constraints, providing a unique and accurate physical index for the subsequent generation of the spatial variation point diffusion function.

[0095] Further, in step S2, using the viewpoint-related optical path depth mapping map generated in step S1 as an index, a preset depth-scattering response model is queried to independently map the corresponding point spread function for each pixel on the imaging plane, and a spatially variable point spread function array is constructed in which the diffusion scale parameter changes nonlinearly with the effective scattering optical path.

[0096] The core of this step is to transform the physical distance (effective scattering optical path) calculated in step S1 into an optical operator that describes the blurring pattern of the image, thereby providing a precise mathematical tool for subsequent image restoration.

[0097] In step S2, the depth-scattering response model uses the logical stearic scattering saturation calculation logic based on physical optics to establish a nonlinear mapping relationship between the effective scattering optical path and the dispersion scale parameter.

[0098] Specifically, traditional Gaussian blur models typically assume a simple linear relationship between the blur radius and depth. This may work when dealing with shallow tissues, but in imaging deep liver tissues, photons undergo multiple scattering and tend to become isotropic after penetrating a certain depth. As a result, the diffusion rate of the blur spot gradually slows down and reaches a physical limit.

[0099] In order to accurately simulate this bio-optical phenomenon, this embodiment abandons the linear model and instead adopts the S-curve growth model, namely the logistic function and its variants, to describe the change of spot size with depth.

[0100] First, we set the system diffraction limit base value, which characterizes the inherent optical properties of the endoscope lens, and the maximum scattering saturation threshold, which characterizes the depletion of photon energy in deep tissues.

[0101] In this step, the system initializes two key boundary parameters.

[0102] Regarding the system's diffraction-limited baseline value, this parameter represents the inherent minimum blur radius caused solely by diffraction from the endoscope lens's optical aperture and sensor pixel sampling when the effective scattered optical path is zero, i.e., the imaging target is located at the outermost surface of the tissue. In this embodiment, this value is preferably set to 0.5 to 1.0 pixel units, and its value is determined based on the modulation transfer function curve of the laparoscopic imaging system and the Rayleigh criterion. Its function is to serve as a baseline for blur calculation, ensuring that even in the non-scattering region, the generated point spread function will not collapse into an infinitely small single point, thereby avoiding high-frequency ringing artifacts generated by the restoration algorithm.

[0103] Regarding the maximum scattering saturation threshold, this parameter represents the upper limit of the maximum diffusion radius that the point spread function can reach in deep tissues after photon energy depletion or the mean free path of scattering reaches statistical equilibrium. In this embodiment, it is preferably set to twenty to thirty pixel units, and its value is based on the measured data of the maximum spot size of near-infrared light in liver tissue slices with a thickness of five to ten millimeters in Monte Carlo photon transport simulation experiments. Its function is to serve as the physical upper limit of fuzzy calculation, preventing the algorithm from generating excessively large convolution kernels when processing extremely deep regions, which could lead to computational crashes or overly smoothed images.

[0104] During the calculation process, the effective scattered optical path output by step S1 is calculated as the numerical difference between the preset ballistic light to scattered light transition optical path. This numerical difference is then used as an input variable and substituted into the exponential decay function. Combined with the preset scattering phase transition rate factor, the nonlinear growth ratio is calculated.

[0105] The system reads the effective scattered optical path corresponding to the current pixel on the imaging plane and introduces the key parameter of the turning optical path from ballistic light to scattered light.

[0106] In this invention, the ballistic-to-scattered light transition path represents the critical depth at which photons, when propagating within tissue, transition from a quasi-ballistic state dominated by straight-line propagation to a random-walk, multi-scattered state dominated by multiple scattering. In this embodiment, it is set to 1.0 to 2.0 millimeters, and its value is determined based on the physical constant of the mean free path of scattering in liver tissue. Its function is to serve as the inflection point of the nonlinear mapping function, controlling the position of the morphological change of the fuzzy growth curve.

[0107] Meanwhile, the system uses the scattering phase transition rate factor to control the steepness of the curve.

[0108] In this invention, the scattering phase transition rate factor characterizes the rate at which the blurriness deteriorates with increasing depth; in this embodiment, it is set to 1.5 to 2.0 per millimeter, and its value is obtained by fitting the optical scattering characteristic curve of isolated liver tissue; its function is to determine the visual gradient of the image transitioning from a clear area to a blurred area.

[0109] In the specific calculation, the system performs a subtraction operation to calculate the optical path difference, multiplies the difference by the rate factor and takes the negative as the power of the natural exponent to calculate the exponential decay term, and then uses the mathematical form of the logistic function to calculate the nonlinear growth ratio between zero and one.

[0110] Perform a multiplication operation, multiply the maximum scattering saturation threshold by the nonlinear growth ratio, and superimpose the product onto the system's diffraction-limited basis value to obtain the initial dispersion scale parameter.

[0111] The system performs a multiplication operation to calculate the product of the maximum scattering saturation threshold and the nonlinear growth ratio, thus obtaining the fuzzy component increased due to the scattering effect.

[0112] Then, an addition operation is performed to add the component to the system's diffraction-limited basis value.

[0113] This calculation process ensures that the diffusion scale parameter is close to the base value in the shallow region and smoothly approaches the saturation threshold in the deep region, which is in line with physical laws.

[0114] Furthermore, the geometric interaction data between the virtual ray output by the reverse ray tracing logic in step S1 and the liver capsule is further invoked to calculate the cosine value of the incident angle and construct anisotropic morphology correction coefficients. The initial diffusion scale parameters are then multiplied to correct the elliptical spot effect caused by oblique incidence.

[0115] It is worth noting that when light shines obliquely into the tissue surface, the circular light spot will be projected as an ellipse.

[0116] The system calculates the cosine of the incident angle of the virtual ray and takes its reciprocal to construct the anisotropic morphology correction coefficient.

[0117] In this invention, the anisotropic morphology correction coefficient characterizes the stretching ratio of the effective light-receiving area caused by the geometric projection angle; in this embodiment, the value range is greater than or equal to one; its value is based on the projection transformation principle in geometric optics; its function is to amplify and correct the diffusion scale of the oblique projection area.

[0118] The system performs a multiplication operation, multiplying the initial diffusion scale parameter by the coefficient to generate a geometrically corrected parameter.

[0119] The final output is the final dispersion scale parameter that can simulate the phase transition process of photons from a quasi-ballistic state to a multiple scattering state.

[0120] The final diffusion scale parameter precisely describes the characteristic width of the spot formed on the imaging plane at the current pixel location after the point light source is scattered by liver tissue with the same effective scattered optical path length as calculated in step S1.

[0121] In step S2, the process of constructing the spatial variation point diffusion function array performs an energy decoupling operation between the ballistic light component and the diffuse light component.

[0122] To further improve the restoration accuracy, this embodiment does not use a single Gaussian kernel, but instead decomposes the optical signal into two parts:

[0123] One part is ballistic light that has not been scattered or has been scattered very little, and retains the high-frequency edge information of the object;

[0124] The other part is diffused light that has undergone multiple random scatterings, forming a low-frequency background fog.

[0125] The ballistic light weights that retain high-frequency edge information and the diffuse light weights that cause background blurring are calculated based on the effective scattering optical path.

[0126] The system calculates energy distribution based on the Beer-Lambert law.

[0127] The ballistic light weight is set to decrease exponentially with the increase of the effective scattered optical path, and the remaining energy after subtracting the ballistic light weight from the total energy is allocated to the diffuse light weight.

[0128] Specifically, the system calculates the ballistic beam weights using the natural constant as the base and the product of the effective scattering path length and the tissue extinction coefficient as the negative exponent.

[0129] In this invention, ballistic light weight represents the proportion of photon energy that still maintains quasi-linear propagation after penetrating tissue; in this embodiment, its range is set to be between zero and one, rapidly approaching zero with increasing depth, and its value is based on the classical biological tissue optical transmission theory; its effect is to determine the extent to which the restoration algorithm retains the sharpness of the original image.

[0130] Subsequently, the system subtracts the ballistic light weight from the numerical value to obtain the diffuse light weight.

[0131] When synthesizing the point spread function kernel of a single pixel, a ballistic light core representing a sharp signal is constructed using a unit impulse function, and a long-tailed distribution smooth function representing a diffuse halo is constructed using a diffusion scale parameter calculated in the above steps.

[0132] The system first constructs an ideal unit impulse function (i.e., a Dirac function matrix with a center pixel of 1 and the rest of the pixels of 0) as the convolution kernel for the ballistic light components to simulate ideal imaging.

[0133] Secondly, the system uses the final dispersion scale parameter obtained from the aforementioned calculation to generate a long-tailed distribution smoothing function.

[0134] In this invention, compared with the Gaussian function, the long-tailed distribution smoothing function (such as the Lorentz function or the Moffat function) can more accurately describe the significant lateral scattering phenomenon in biological tissues. In this embodiment, the long-tailed distribution smoothing function is taken as its full width at half maximum (FWHM) is strictly equal to the final diffusion scale parameter, and its value is based on the experimentally measured waveform fitting results of the point spread function. Its effect is to accurately simulate the diffusion morphology of the halo.

[0135] Then, the two are multiplied by their respective weights and then linearly superimposed.

[0136] The system performs matrix addition and scalar multiplication operations, multiplying the ballistic light core by the ballistic light weight, multiplying the diffuse light periphery by the diffuse light weight, and adding the two together to synthesize the unique hybrid point diffusion function kernel of the current pixel.

[0137] Traverse all pixel positions on the imaging plane, and spatially tile and tensor reassemble each generated point spread function kernel with independent size and shape according to pixel coordinates.

[0138] The system repeats the above calculation process for each coordinate point within the image resolution range.

[0139] The final result is a spatial variation point spread function array with a four-dimensional structure. Each element in the spatial variation point spread function array corresponds precisely to the optical degradation characteristics of that location in the image after being affected by physical scattering.

[0140] The array is a four-dimensional data volume containing height, width, kernel height, and kernel width. It provides a point-by-point accurate physical operator for the inverse deconvolution in step S3, ensuring that the image restoration process is based on the inverse operation of the real physical process, rather than blind image enhancement.

[0141] Further, in step S3, the spatial variation point diffusion function array constructed in step S2 is called as the inverse restoration operator, and the central skeleton projection line of the MRCP three-dimensional biliary tract model in step S1 is used as the topological gravity constraint to perform restricted inverse deconvolution processing on the two-dimensional ICG fluorescence image obtained in step S1, thereby eliminating scattered halo and generating a restored biliary tract image.

[0142] This step utilizes the physical operator (point spread function array) generated in step S2 and the anatomical reference (MRCP model) generated in step S1 to restore the blurred two-dimensional ICG fluorescence image into a clear and accurate image of the biliary tract structure.

[0143] In step S3, the specific logic for establishing topological gravitational constraints employs the technique of constructing a dissected potential energy field.

[0144] Specifically, the system uses geometric data stored in computer memory to construct a two-dimensional scalar field with a resolution that is completely consistent with that of the imaging plane.

[0145] Each value in this scalar field no longer represents brightness, but rather the prior probability that a real biliary signal exists at that location. In this way, the system transforms the hard geometric constraints of the MRCP 3D model into a soft mathematical potential field, thus enabling smooth integration into subsequent iterative restoration algorithms.

[0146] First, extract the central skeleton projection line generated by projecting the MRCP three-dimensional bile duct model through reverse ray tracing logic in step S1, and use it as the standard spatial reference path for anatomy.

[0147] During implementation, the system processor reads the central skeleton data from the output buffer of step S1. This data is usually composed of a series of ordered two-dimensional floating-point coordinates or parameterized spline curves. It accurately depicts the theoretical projection center axis of the common hepatic duct, common bile duct and left and right hepatic ducts on the imaging sensor plane from the current laparoscopic perspective. This path is defined as the theoretical geometric true value of the anatomical structure and serves as the absolute spatial reference benchmark for all subsequent topological constraint calculations.

[0148] Then, for each pixel in the resolution grid of the imaging plane, the Euclidean distance from that pixel to the nearest central skeleton projection line segment is calculated, and this Euclidean distance is defined as the skeleton Euclidean deviation.

[0149] The system starts a parallel computing thread to traverse the coordinates of each pixel in the image matrix. For the currently processed pixel, the system performs the minimum Euclidean distance calculation based on the principle of analytical geometry. That is, it calculates the vertical projection distance or endpoint distance from the pixel coordinates to each line segment on the central skeleton projection line, and selects the minimum value as the measurement result of the point. This calculation result is stored in a distance mapping matrix of the same size in real time.

[0150] The Euclidean deviation of the skeleton characterizes the geometric distance of the current fluorescent signal point from the ideal anatomical path.

[0151] From a physical perspective, if the value is zero, it means that the pixel is precisely located at the theoretical center of the bile duct. If the value increases monotonically as the pixel moves away from the center line, it means that the probability that the position belongs to background noise or non-specific signal gradually increases.

[0152] A potential energy attenuation weight field is constructed based on a preset anatomical tolerance radius. This anatomical tolerance radius defines the maximum allowable spatial registration error limit between the MRCP three-dimensional biliary tract model and the real biliary tract in the surgical scenario.

[0153] In this step, the system introduces a crucial spatial threshold parameter: the anatomical tolerance radius.

[0154] The anatomical tolerance radius represents the comprehensive confidence interval radius of the system for non-rigid deformation of soft tissue, organ displacement caused by respiratory motion, and initial registration residuals. Geometrically, it defines the boundary radius threshold of the high-confidence spatial tolerance region. In this embodiment, the radius is preferably set to 3.0 millimeters. During specific calculations, the system converts this physical length into the corresponding pixel value according to the imaging resolution of the camera (the number of pixels per millimeter). For example, it is about 30 to 50 pixels under a high-definition endoscope. Its value is derived from the clinical accuracy standard of hepatobiliary surgical navigation systems. Statistical data show that after non-rigid registration, an error within 3 millimeters is within the acceptable range for surgical safety, while deviations exceeding this range usually indicate registration failure or signal abnormalities.

[0155] The anatomical tolerance radius determines the geometric width of the maximum weighted response interval of the potential energy field. Within this radius, the algorithm gives all signals a very high degree of confidence, allowing them to be enhanced and restored; ensuring that even if real biliary signals undergo slight displacement, they will not be misjudged as noise and rejected.

[0156] When calculating the potential energy weight, the ratio of the Euclidean deviation of the skeleton to the anatomical tolerance radius is used as input and substituted into a high-order reciprocal polynomial function with the gravitational field decay order as the power, to calculate the potential energy decay weight value between zero and one.

[0157] During the calculation, a division operation is first performed to calculate the quotient of the skeleton Euclidean deviation of the current pixel and the anatomical tolerance radius to obtain the normalized distance. Then, the shape control parameter of the gravitational field attenuation order is introduced.

[0158] The gravitational field attenuation order is a dimensionless positive even integer used to control the steepness of the potential field edge roll-off characteristics. In this embodiment, it is preferably set to six, which is based on the design principle of Butterworth filters in signal processing. Compared with second-order attenuation or Gaussian attenuation, sixth-order attenuation can provide a flatter passband and a steeper stopband.

[0159] The high order of the gravitational field decay order ensures that within the tolerance radius, the weight value can be maintained at a level close to one (e.g., 0.99) for a long time. Once it exceeds the radius, the weight value decays to zero at an extremely fast rate (sixth power). This avoids premature suppression of effective signals at the edge and also achieves a decisive cut-off of far-end noise.

[0160] Perform an exponentiation operation to calculate the sixth power of the normalized distance; then add the exponent value to the numerical value; finally, calculate the reciprocal of the sum to obtain the final potential energy decay weight value.

[0161] This computational logic enables the potential energy decay weight field to reach its peak on the central skeleton projection line and maintain a flat high-weight plateau within the anatomical tolerance radius. Beyond this radius, it rapidly decays to zero, thus forming a high-weight distribution area along the skeleton line in geometric space. This defines a high-confidence spatial range for the subsequent deconvolution iteration process and severely suppresses non-specific background fluorescence signals or vascular leakage interference located outside this range.

[0162] Through the above steps, the system mathematically constructs a probability constraint distribution structure with high spatial selectivity, specifying the effective region of signal convergence based on anatomical priors for subsequent iterative algorithms.

[0163] In step S3, the restricted inverse deconvolution process employs a dual-constraint iterative update mechanism based on tensor product.

[0164] This embodiment employs an improved spatial variation Richardson-Lucy iterative algorithm architecture, which extends the traditional one-dimensional convolution into a tensor product operation adapted to the heterogeneity of biological tissues.

[0165] In each iteration cycle, the spatial variation point spread function array constructed in step S2 is first called to perform forward physical scattering simulation on the currently estimated restored image. That is, for each pixel position, its own point spread function kernel is used for local integration to generate a simulated scattering image.

[0166] The system performs spatial mutation convolution. Since the blurring degree varies at different locations in the image, the traditional fast Fourier transform convolution is no longer applicable. The system traverses each pixel in the image, uses the two-dimensional coordinates of the pixel as an index, and addresses and extracts the point spread function kernel that has a unique mapping relationship with the coordinate position from the four-dimensional array generated in step S2.

[0167] Subsequently, with the pixel as the center, the kernel obtained based on coordinate index is locally weighted and summed with the estimated image of the current iteration step. This process mathematically performs a forward degenerate projection operation based on a physical model, that is, using the estimated image of the current iteration step as a potential signal source, the theoretical observed brightness distribution on the imaging plane after being scattered and modulated by non-uniform medium is calculated.

[0168] Subsequently, the brightness ratio of the two-dimensional ICG fluorescence image obtained in step S1 to the simulated scattering image is calculated, and this brightness ratio is back-projected back into the source space through the accompanying inverse operator to construct the physical likelihood term.

[0169] The system performs pixel-by-pixel division, dividing the actually observed two-dimensional ICG fluorescence image by the generated simulated scattering image to obtain the error ratio matrix.

[0170] If the simulation is completely accurate, the ratio should be one everywhere. Then, the transpose matrix of the point spread function array (i.e., the adjoint operator) is constructed, and the transpose matrix is ​​used to perform a reverse convolution operation on the error ratio matrix. The purpose of this step is to back-project the residual between the observed data and the physical model back to the source image space according to the principle of optical path reversibility, and calculate the physical likelihood gain factor used to correct the estimated image.

[0171] Based on this, gradient direction consistency constraint logic is introduced to calculate the brightness gradient direction vector of the restored image in the current iteration step, and to calculate the cosine similarity of the angle between the brightness gradient direction vector and the normal vector of the central skeleton projection line, thus constructing the gradient skeleton consistency tensor.

[0172] To further eliminate noise that is within the anatomical range but does not conform to the shape (such as clump noise), geometric morphological constraints are introduced, and the gradient vector field of the current restored image is calculated using the Sobel operator.

[0173] Subsequently, the dot product (i.e., cosine similarity) of the gradient vector and the normal vector of the MRCP skeleton projection line is calculated.

[0174] The generated gradient skeleton consistency tensor is a normalized scalar field that characterizes the degree of geometric consistency between the direction of the brightness gradient change of the current fluorescence signal and the direction of the theoretical skeleton normal. In this embodiment, its value is between zero and one. Its value is based on the differential geometry of tubular objects. In the real biliary structure, the direction of its brightness gradient should be strictly perpendicular to the tube wall, that is, parallel to the skeleton normal. Therefore, when the gradient direction is consistent with the skeleton normal, the value is one, and the algorithm will enhance the signal; when the gradient direction is disordered (such as Gaussian white noise), the value approaches zero, and the algorithm will suppress the signal.

[0175] Finally, the physical likelihood term is combined with the potential energy decay weight field and gradient skeleton consistency tensor from the above steps by performing a joint multiplication operation, and a topological penalty factor is introduced to adjust the intensity of the intervention of the anatomical prior on the physical reconstruction. Through multiple rounds of iterative updates, the fluorescence signal energy is forced to move towards the skeleton in spatial position and conform to the tubular structure characteristics in morphology until the relative entropy change rate between two adjacent iterations is lower than the preset threshold.

[0176] Specifically, the estimated image from the previous round is multiplied by three factors in sequence: the physical likelihood term, the potential decay weight field, and the gradient skeleton consistency tensor.

[0177] The latter two terms are constraint terms, which need to be multiplied first and then exponentialized. Their exponents are the topological penalty factors.

[0178] The topology penalty factor is a weighted hyperparameter that adjusts the anatomical model constraints relative to the fidelity of the physical data. In this embodiment, it is set to 1.0 to 1.5. Its value is based on an empirical value of adaptive adjustment based on the image signal-to-noise ratio. The lower the signal-to-noise ratio, the larger the value of this factor. The role of the topology penalty factor is that the larger the topology penalty factor, the closer the restoration result is to the shape of the MRCP model; the smaller the topology penalty factor, the more the restoration result depends on the original fluorescence data.

[0179] The above steps are repeated, and the relative entropy (Kullback-Leibler divergence) of the results of two adjacent iterations is calculated at the end of each round. When the rate of change of the relative entropy is lower than a preset threshold (e.g., one ten-thousandth), the system determines that the algorithm has converged, stops the iteration, and outputs the final high-definition restored bile duct image.

[0180] Further, in step S4, the anatomical contour features of the reconstructed bile duct image generated in step S3 are extracted, the spatial registration relationship between the MRCP three-dimensional bile duct model and the reconstructed bile duct image is corrected according to the anatomical contour features, and the MRCP three-dimensional bile duct model is rendered and superimposed onto the reconstructed bile duct image.

[0181] Step S4 constructs a closed-loop registration correction system, using the high-quality image restored in step S3 as the observation ground truth, to perform non-rigid calibration on the projection of the MRCP model in screen space, and establishes a confidence-based secure display mechanism.

[0182] In step S4, the specific logic for performing geometric consistency verification includes the construction process of the local confidence graph.

[0183] To prevent residual artifacts in the restored image from misleading the registration algorithm, it is first necessary to quantify the reliability of image features.

[0184] First, a multi-scale Hessian matrix filter is applied to the reconstructed bile duct image generated in step S3 to extract response feature maps that can enhance the contrast of tubular structures and suppress sheet noise. Then, a reconstructed skeleton binary map reflecting the centerline of the tubular structure is generated through binarization.

[0185] Multi-scale Hessian matrix analysis is employed. The eigenvalues ​​of the Hessian matrix reflect the local second-order derivative structure of the image, i.e., curvature information. In this embodiment, the scale parameter is set to one to four, covering bile ducts of different diameters, and its value is based on the pixel mapping range of the anatomical diameter of the bile duct. This method specifically enhances tubular structures (linear structures) while suppressing patchy reflections (spotted structures) on the liver surface.

[0186] Subsequently, an adaptive threshold segmentation algorithm is applied to transform the response feature map into a single-pixel-width restored skeleton binary map, where the pixel value on the skeleton is one and the background value is zero.

[0187] Simultaneously, based on the real-time pose matrix obtained in step S1, a binary projection skeleton map of the MRCP 3D biliary tract model is generated at the current viewpoint.

[0188] Using the current camera extrinsic matrix, the centerline data of the MRCP model is projected onto the imaging plane to generate a binary projection skeleton map. These two binary maps respectively represent the measured anatomical topology based on intraoperative real-time fluorescence signal reconstruction and the theoretical anatomical topology based on the rigid pose projection of the preoperative model.

[0189] Subsequently, a local sliding window centered on the current pixel is established, and the normalized cross-correlation coefficient between the restored skeleton binary map and the projected skeleton binary map is calculated within the local sliding window to quantify the local similarity between the two in terms of topological morphology.

[0190] The local sliding window defines the local receptive field for texture matching. In this embodiment, it is set to 16 x 16 pixels, and its value is based on empirical values. This size is sufficient to include the bifurcation features of the bile duct. Its function is to calculate the statistical correlation between two binary images within the local sliding window. If the two structures overlap, the normalized cross-correlation coefficient approaches one; if they are misaligned or do not match in shape, the coefficient approaches zero or a negative value.

[0191] Based on this, the gradient magnitude of the reconstructed biliary tract image is further calculated and a gradient significance factor is constructed. Multiplication is performed, and the gradient significance factor is jointly weighted with the normalized cross-correlation coefficient to generate a local confidence map.

[0192] The gradient significance factor is used to suppress the modulation function of the weights in flat regions (featureless regions). In this embodiment, the hyperbolic tangent function is used for nonlinear mapping. The logic is to calculate the gradient magnitude divided by the noise base threshold (e.g., 5.0) and substitute the quotient into the hyperbolic tangent function. Its value is based on the fact that the hyperbolic tangent function has good normalization properties, which can map any positive number to between zero and one, and has a linear response near the zero point. Its function is that when the gradient magnitude is lower than the noise base, the factor output is close to zero; when the gradient magnitude is significant, the factor smoothly saturates to one.

[0193] Then, a multiplication operation is performed to multiply the gradient significance factor by the normalized cross-correlation coefficient to obtain the final local confidence plot.

[0194] Each value in this local confidence map accurately represents the degree of topological fit between the image features and model features at the corresponding pixel location, and the confidence of flat regions with low image gradient magnitudes is forced to zero by introducing a gradient significance factor.

[0195] This step, through gradient weighting logic, mathematically suppresses regions with high cross-correlation responses but insufficient gradient magnitude significance that lack geometric constraints, thereby preventing the iterative nearest point algorithm from generating incorrect point set correspondences or numerical convergence deviations in low-contrast homogeneous regions due to a lack of feature constraints.

[0196] In step S4, the logic for correcting spatial registration relationships and rendering overlay adopts a view plane non-rigid deformation compensation and visual blocking rendering mechanism.

[0197] This embodiment uses screen space distortion technology to correct projection errors.

[0198] First, in the local confidence map generated in the above steps, high confidence regions with values ​​higher than the preset safety threshold are selected. Then, in the high confidence region, the iterative nearest neighbor algorithm is used to search for the nearest neighbor correspondence between the skeleton point set of the restored biliary tract image and the projection skeleton point set of the MRCP three-dimensional biliary tract model. The planar displacement difference between each pair of feature points is calculated to construct the residual deviation vector field.

[0199] In this embodiment, the preset safety threshold is set to 0.75, which serves as a high-pass filter for feature point screening.

[0200] Within the selected region, point pairs are established using the iterative nearest-point algorithm, and the residual deviation vector is calculated. This residual deviation vector is a two-dimensional vector that indicates the direction and distance that the model projection points need to move to coincide with the restored image.

[0201] Subsequently, the residual deviation vector field is fitted globally using a thin plate spline transformation model to construct a two-dimensional mapping function that can describe the nonlinear deformation of the entire field. This two-dimensional mapping function is then used to update the screen output coordinates of the mesh vertices of the MRCP three-dimensional biliary tract model in the rendering pipeline, thereby achieving non-rigid deformation compensation only for the view plane projection shape without modifying the physical coordinates of the three-dimensional model.

[0202] Thin Plate Spline (TPS) transformation function is an interpolation method based on radial basis functions to simulate the bending shape of an infinitely large thin metal plate under the action of constraint points. In this embodiment, the kernel function is the square of the distance multiplied by the natural logarithm of the distance. The value is derived from the fact that thin plate splines can ensure the continuity of the second derivative of the deformation field, that is, the bending energy is minimized, which is highly consistent with the smooth deformation characteristics of biological soft tissues. Its function is to generalize the discrete deviation vector into a continuous deformation field across the entire screen.

[0203] In the graphics rendering pipeline, the system writes vertex shader programs.

[0204] The program receives the original vertices of the MRCP model, projects them onto screen coordinates, adds the local offset calculated by the TPS function, and outputs the corrected screen coordinates.

[0205] This method avoids modifying the physical coordinates of 3D vertices, thus ensuring the continuity of the depth buffer.

[0206] In the final rendering stage, the system introduces hierarchical display logic. For areas with a confidence level higher than the safety threshold, the model outline after deformation compensation is drawn with a highlighted solid line.

[0207] For high-confidence areas, the system outputs pure green pixels (RGB value 0,255,0) with a line width set to two to three pixels to provide clear navigation.

[0208] For areas with confidence levels below the safety threshold, the system performs pixel discarding to prevent rendering of the model outline and triggers a visual blocking mechanism, generating a dynamically flashing semi-transparent warning cloud map to cover the area, so as to intuitively remind the operator that there is registration uncertainty in this area.

[0209] For low-confidence regions, the system executes a fragment discard instruction and does not draw the model outline. Instead, a semi-transparent warning cloud is rendered in the post-processing stage.

[0210] The specific implementation logic of dynamic blinking is as follows: The transparency channel of the pixel is controlled by a time modulation function. In this embodiment, the transparency value is equal to the base transparency of 0.3 plus the modulation amplitude of 0.2 multiplied by the sine function. The independent variable of the sine function is the system time multiplied by the frequency factor (such as three radians per second). Its function is to make the red warning area present a periodic light and dark breathing effect (i.e., blinking). By taking advantage of the human visual system's sensitivity to dynamic changes, it forces the doctor's attention and warns him that the area is unreliable.

[0211] Example 2

[0212] like Figure 2 As shown, the present invention also discloses an ICG fluorescence image and MRCP three-dimensional model registration and fusion display system, comprising: a viewpoint-related optical path depth map generation module, a spatial variation point diffusion function array construction module, a restricted inverse deconvolution restoration module, and a registration and fusion display module, wherein;

[0213] The viewpoint-related optical path depth mapping generation module is configured to simultaneously acquire the preoperatively reconstructed MRCP three-dimensional biliary tract model, the intraoperatively acquired two-dimensional ICG fluorescence image and the corresponding laparoscopic camera pose matrix, perform viewpoint projection on the MRCP three-dimensional biliary tract model based on the pose matrix, and calculate the effective scattering optical path of the anatomical structure corresponding to each pixel on the imaging plane in the liver tissue through reverse ray tracing logic, and generate a viewpoint-related optical path depth mapping containing tissue thickness information.

[0214] The spatial variation point diffusion function array construction module is configured to use the viewpoint-related optical path depth mapping generated by the viewpoint-related optical path depth mapping generation module as an index to query the preset depth-scattering response model, independently map the corresponding point diffusion function for each pixel on the imaging plane, and construct a spatial variation point diffusion function array whose diffusion scale parameter changes nonlinearly with the effective scattering optical path.

[0215] The constrained inverse deconvolution restoration module is configured to call the spatial variation point diffusion function array constructed by the spatial variation point diffusion function array construction module as the inverse restoration operator, and use the central skeleton projection line of the MRCP three-dimensional biliary tract model generated by the viewpoint-related optical path depth map generation module as the topological gravity constraint to perform constrained inverse deconvolution processing on the two-dimensional ICG fluorescence image obtained by the viewpoint-related optical path depth map generation module to eliminate scattered halo and generate restored biliary tract image;

[0216] The registration and fusion display module is configured to extract the anatomical contour features of the restored bile duct image generated by the restricted inverse deconvolution restoration module, correct the spatial registration relationship between the MRCP three-dimensional bile duct model and the restored bile duct image based on the anatomical contour features, and render and overlay the MRCP three-dimensional bile duct model onto the restored bile duct image.

[0217] The above description of the disclosed embodiments enables those skilled in the art to make or use the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for registering and fusing ICG fluorescence images with MRCP three-dimensional models, characterized in that, Includes the following steps: Step S1: Simultaneously acquire the preoperative reconstructed MRCP three-dimensional biliary tract model, the intraoperative real-time acquired two-dimensional ICG fluorescence image and the corresponding laparoscopic camera pose matrix. Based on the pose matrix, perform viewpoint projection on the MRCP three-dimensional biliary tract model. Calculate the effective scattering optical path of the anatomical structure corresponding to each pixel on the imaging plane in the liver tissue through reverse ray tracing logic, and generate a viewpoint-related optical path depth mapping map containing tissue thickness information. Step S2: Using the viewpoint-related optical path depth mapping map generated in step S1 as an index, query the preset depth-scattering response model, independently map the corresponding point spread function for each pixel on the imaging plane, and construct a spatially variable point spread function array in which the diffusion scale parameter changes nonlinearly with the effective scattering optical path. Step S3: Call the spatial variation point diffusion function array constructed in step S2 as the inverse restoration operator, and use the central skeleton projection line of the MRCP three-dimensional biliary tract model in step S1 as the topological gravity constraint to perform restricted inverse deconvolution processing on the two-dimensional ICG fluorescence image obtained in step S1, eliminate the scattered halo and generate the restored biliary tract image. Step S4: Extract the anatomical contour features of the reconstructed bile duct image generated in step S3, correct the spatial registration relationship between the MRCP three-dimensional bile duct model and the reconstructed bile duct image based on the anatomical contour features, and render and overlay the MRCP three-dimensional bile duct model onto the reconstructed bile duct image.

2. The method for registering, fusing, and displaying ICG fluorescence images and MRCP three-dimensional models according to claim 1, characterized in that, In step S1, the specific execution process of the reverse ray tracing logic includes: Using the optical center of the laparoscopic camera as the origin, virtual rays are emitted to each pixel on the imaging plane to perform multi-level collision detection on the MRCP three-dimensional biliary model. Calculate the first intersection point between the virtual ray and the liver capsule mesh in the MRCP three-dimensional biliary tract model and mark it as the media incident point; Calculate the second intersection point of the same ray with the intrahepatic bile duct grid and mark it as the scattering termination point; Offset sampling is performed along the path connecting the incident point and the scattering termination point of the medium. Based on the pre-annotated anatomical semantic labels in the MRCP three-dimensional biliary model, non-scattering medium regions belonging to portal vein vessels or liver cysts are identified along the path. The length of the non-scattering medium region is subtracted from the Euclidean distance between the incident point and the scattering termination point to obtain the pure geometric penetration depth containing only highly scattered liver parenchyma components.

3. The method for registering, fusing, and displaying ICG fluorescence images and MRCP three-dimensional models according to claim 2, characterized in that, In step S1, during the generation of the viewpoint-related optical path depth map, the calculation of the effective scattered optical path includes weighted correction logic based on patient physiological heterogeneity: Obtain the patient's preoperative magnetic resonance imaging data or body mass index, extract the patient-specific fat scattering coefficient that characterizes liver fat content, and the normalized index of magnetic resonance signal intensity that characterizes the density of local liver tissue. Physiological scattering enhancement factors were constructed using patient-specific fat scattering coefficients and normalized magnetic resonance signal intensity indices. These factors were then multiplicatively weighted to the obtained pure geometric penetration depth to reflect the different effects of different fat densities on photon scattering behavior. The angle between the virtual ray and the normal vector of the liver capsule mesh at the incident point of the medium is calculated. A Fresnel incident efficiency constraint function is introduced to numerically compensate for the light energy loss caused by large-angle grazing. Finally, the value after geometric semantic elimination, physiological density weighting and incident angle correction is determined as the effective scattered optical path.

4. The method for registering, fusing, and displaying ICG fluorescence images and MRCP three-dimensional models according to claim 3, characterized in that, In step S2, the depth-scattering response model employs a logistic scattering saturation calculation logic based on physical optics to establish a nonlinear mapping relationship between the effective scattering optical path and the dispersion scale parameter: Set the system diffraction limit base value to characterize the inherent optical properties of the endoscope lens, and the maximum scattering saturation threshold to characterize the photon energy depletion of deep tissue; The numerical difference between the effective scattered optical path output in step S1 and the preset ballistic light to scattered light turning optical path is calculated, and this numerical difference is used as an input variable and substituted into the exponential decay function. The nonlinear growth ratio is calculated in combination with the preset scattering phase transition rate factor. The initial dispersion scale parameter is obtained by multiplying the maximum scattering saturation threshold by the nonlinear growth ratio and then superimposing the product onto the system diffraction-limited basis value. The geometric interaction data between the virtual ray output by the reverse ray tracing logic in step S1 and the liver capsule is called to calculate the cosine value of the incident angle and construct the anisotropic morphology correction coefficient. The initial diffusion scale parameter is multiplied to correct the elliptical spot effect caused by oblique incidence. The output is the final dispersion scale parameter that can simulate the phase transition process of photons from a quasi-ballistic state to a multiple scattering state.

5. The method for registering, fusing, and displaying ICG fluorescence images and MRCP three-dimensional models according to claim 4, characterized in that, In step S2, the process of constructing the spatial variation point diffusion function array performs an energy decoupling operation between the ballistic light component and the diffuse light component: The ballistic beam weights that retain high-frequency edge information and the diffuse beam weights that cause background blur are calculated based on the effective scattered optical path. The ballistic light weight is set to decrease exponentially with the increase of the effective scattered optical path, and the remaining energy after subtracting the ballistic light weight from the total energy is allocated to the diffuse light weight. When synthesizing the point spread function kernel of a single pixel, the ballistic light core representing the sharp signal is constructed using the unit impulse function, and the long-tailed distribution smooth function generated by the calculated dispersion scale parameter is used to construct the diffuse light periphery representing the scattered halo. Then, the two are multiplied by their respective weights and then linearly superimposed. Traverse all pixel positions on the imaging plane, and spatially tile and tensor reassemble each generated point spread function kernel with independent size and shape according to pixel coordinates; The final result is a spatial variation point diffusion function array with a four-dimensional structure. Each element in the spatial variation point diffusion function array corresponds precisely to the optical degradation characteristics of that location in the image after being affected by physical scattering.

6. The method for registering, fusing, and displaying ICG fluorescence images and MRCP three-dimensional models according to claim 5, characterized in that, In step S3, the topological gravitational constraint is established using the dissected potential field construction technique: Extract the central skeleton projection line generated by projecting the MRCP three-dimensional bile duct model through reverse ray tracing logic in step S1, and use it as the standard spatial reference path for anatomy. For each pixel in the imaging plane resolution grid, the Euclidean distance from that pixel to the nearest central skeleton projection line segment is calculated, and this Euclidean distance is defined as the skeleton Euclidean deviation; the skeleton Euclidean deviation characterizes the geometric distance of the current fluorescence signal point from the ideal anatomical path; A potential energy attenuation weight field is constructed based on the preset anatomical tolerance radius. This anatomical tolerance radius defines the maximum allowable spatial registration error limit between the MRCP three-dimensional biliary tract model and the real biliary tract in the surgical scenario. The ratio of the Euclidean deviation of the skeleton to the anatomical tolerance radius is used as input and substituted into a high-order reciprocal polynomial function with the gravitational field decay order as the power to calculate the potential energy decay weight value between zero and one.

7. The method for registering, fusing, and displaying ICG fluorescence images and MRCP three-dimensional models according to claim 6, characterized in that, In step S3, the constrained inverse deconvolution process employs a dual-constraint iterative update mechanism based on tensor product: In each iteration cycle, the spatial variation point spread function array constructed in step S2 is first called to perform forward physical scattering simulation on the currently estimated restored image. That is, for each pixel position, its own point spread function kernel is used for local integration to generate a simulated scattering image. The brightness ratio of the two-dimensional ICG fluorescence image and the simulated scattering image obtained in step S1 is calculated, and the brightness ratio is back-projected back to the source space through the adjoint inverse operator to construct the physical likelihood term; Introducing gradient direction consistency constraint logic, the brightness gradient direction vector of the restored image in the current iteration step is calculated, and the cosine similarity of the angle between the brightness gradient direction vector and the normal vector of the central skeleton projection line is calculated to construct the gradient skeleton consistency tensor. The physical likelihood term, the potential energy decay weight field, and the gradient skeleton consistency tensor are subjected to joint multiplication. A topological penalty factor is introduced to adjust the intensity of the intervention of the anatomical prior on the physical reconstruction. Through multiple rounds of iterative updates, the fluorescence signal energy is forced to move towards the skeleton in spatial position and conform to the tubular structure characteristics in morphology until the relative entropy change rate between two adjacent iterations is lower than a preset threshold.

8. The method for registering, fusing, and displaying ICG fluorescence images and MRCP three-dimensional models according to claim 7, characterized in that, In step S4, the specific logic for performing geometric consistency verification includes the construction process of the local confidence graph: A multi-scale Hessian matrix filter is applied to the reconstructed bile duct image generated in step S3 to extract response feature maps that can enhance the contrast of tubular structures and suppress sheet noise. The reconstructed skeleton binary map reflecting the center line of the tubular structure is generated through binarization. Based on the real-time pose matrix obtained in step S1, a binary projection skeleton map of the MRCP 3D biliary tract model is generated at the current viewpoint. Establish a local sliding window centered on the current pixel, and calculate the normalized cross-correlation coefficient between the restored skeleton binary map and the projected skeleton binary map within the local sliding window to quantify the local similarity of the two in terms of topological morphology. The gradient magnitude of the reconstructed bile duct image is calculated and a gradient significance factor is constructed. Multiplication is performed, and the gradient significance factor is jointly weighted with the normalized cross-correlation coefficient to generate a local confidence map. Each value in the local confidence map accurately represents the degree of topological fit between the image features and model features at the corresponding pixel location. Furthermore, by introducing the gradient significance factor, the confidence of flat regions with excessively low image gradient magnitudes is forced to be reduced to zero.

9. The method for registering, fusing, and displaying ICG fluorescence images and MRCP three-dimensional models according to claim 8, characterized in that, In step S4, the logic for correcting spatial registration and rendering overlay employs a view plane non-rigid deformation compensation and visual blocking rendering mechanism: In the local confidence map, high confidence regions with values ​​higher than the preset safety threshold are selected, and the iterative nearest point algorithm is used in the high confidence region to search for the nearest neighbor correspondence between the skeleton point set of the restored biliary tract image and the projection skeleton point set of the MRCP three-dimensional biliary tract model. The planar displacement difference between each pair of feature points is calculated to construct the residual deviation vector field. The residual deviation vector field is fitted globally using a thin plate spline transformation model to construct a two-dimensional mapping function that can describe the nonlinear deformation of the entire field. This two-dimensional mapping function is then used to update the screen output coordinates of the mesh vertices of the MRCP three-dimensional biliary tract model in the rendering pipeline. In the final rendering stage, a hierarchical display logic is introduced. For areas with a confidence level higher than the safety threshold, the model outline after deformation compensation is drawn with a highlighted solid line. For areas with confidence levels below the safety threshold, a pixel discard operation is performed to prevent the rendering of the model outline, and a visual blocking mechanism is triggered to generate a dynamically flashing semi-transparent warning cloud map to cover the area, so as to intuitively remind the operator that there is registration uncertainty in this area.

10. A system for registering and fusing ICG fluorescence images with MRCP three-dimensional models, based on the method for registering and fusing ICG fluorescence images with MRCP three-dimensional models according to any one of claims 1-9, characterized in that, include: The system includes a viewpoint-related optical path depth mapping generation module, a spatial variation point spread function array construction module, a restricted inverse deconvolution restoration module, and a registration and fusion display module. The viewpoint-related optical path depth mapping generation module is configured to simultaneously acquire the preoperatively reconstructed MRCP three-dimensional biliary tract model, the intraoperatively acquired two-dimensional ICG fluorescence image and the corresponding laparoscopic camera pose matrix, perform viewpoint projection on the MRCP three-dimensional biliary tract model based on the pose matrix, and calculate the effective scattering optical path of the anatomical structure corresponding to each pixel on the imaging plane in the liver tissue through reverse ray tracing logic, and generate a viewpoint-related optical path depth mapping containing tissue thickness information. The spatial variation point diffusion function array construction module is configured to use the viewpoint-related optical path depth mapping generated by the viewpoint-related optical path depth mapping generation module as an index to query the preset depth-scattering response model, independently map the corresponding point diffusion function for each pixel on the imaging plane, and construct a spatial variation point diffusion function array whose diffusion scale parameter changes nonlinearly with the effective scattering optical path. The constrained inverse deconvolution restoration module is configured to call the spatial variation point diffusion function array constructed by the spatial variation point diffusion function array construction module as the inverse restoration operator, and use the central skeleton projection line of the MRCP three-dimensional biliary tract model generated by the viewpoint-related optical path depth map generation module as the topological gravity constraint to perform constrained inverse deconvolution processing on the two-dimensional ICG fluorescence image obtained by the viewpoint-related optical path depth map generation module to eliminate scattered halo and generate restored biliary tract image; The registration and fusion display module is configured to extract the anatomical contour features of the restored bile duct image generated by the restricted inverse deconvolution restoration module, correct the spatial registration relationship between the MRCP three-dimensional bile duct model and the restored bile duct image based on the anatomical contour features, and render and overlay the MRCP three-dimensional bile duct model onto the restored bile duct image.

Images