A time-domain diffuse optical tomography reconstruction method and system

Abstract

The application belongs to the field of near-infrared spectroscopy and biomedical engineering, and specifically discloses a time-domain diffuse optical tomography reconstruction method and system. The method is as follows: a three-dimensional geometric model is constructed, and a time-domain diffusion equation is used to describe the transmission process of light in biological tissues; the finite element method is used to solve the equation and obtain the time point diffusion function of the surface of the three-dimensional geometric model; feature data is extracted according to the time point diffusion function, and spatial sensitivity distribution maps of the feature data are drawn respectively; based on the spatial sensitivity distribution maps and node depth, feature data with high spatial sensitivity are selected as the basis for reconstruction; an inverse problem optimization model is established, and the optical parameter reconstruction of shallow nodes and deep nodes is sequentially completed by using the corresponding feature data. The application improves the stability and reliability of optical reconstruction.

Description

Technical Field

[0001] This application belongs to the fields of near-infrared spectroscopy and biomedical engineering, and more specifically, relates to a time-domain diffuse optical tomography reconstruction method and system. Background Technology

[0002] Diffused optical computed tomography (DOT) is an emerging non-invasive imaging technique that can reconstruct optical parameters such as absorption and scattering coefficients within biological tissues. This technique utilizes the low absorption rate of biological tissues in the near-infrared band (600nm~1100nm), allowing for detection depths of several centimeters. Due to its non-invasive and low-cost advantages, it has significant application value in areas such as early diagnosis of breast cancer, neonatal brain functional imaging, and small animal imaging.

[0003] Depending on the type of light source, diffuse optical computed tomography (DCT) can be divided into three modes: continuous wave, frequency domain, and time domain. Among them, the continuous wave system uses continuous laser illumination to measure the light intensity distribution on the tissue surface, and its system structure is simple and the lowest cost; the frequency domain system uses a radio frequency modulated light source to measure the amplitude attenuation and phase delay of the modulated light, and the cost of the light source and detector is moderate; the time domain method uses ultrashort pulse lasers to measure the time-of-flight distribution of photons in the tissue, and although this system is more expensive, it can obtain the richest optical information.

[0004] However, the inverse problem of diffuse optical tomography is not unique. Continuous wave systems only measure light intensity information, lacking time or phase resolution capabilities and having limited optical parameter reconstruction capabilities; frequency domain systems can provide richer measurement information to a certain extent through phase measurement; time domain methods can obtain the most comprehensive information dimensions, thus becoming the most promising technical route for diffuse optical tomography. Making full use of time domain data while minimizing computational complexity has become a key scientific problem in this field that is both challenging and cutting-edge.

[0005] Commonly used data types in time-domain diffuse optical tomography (DOT) include: full-time-resolution data, time-window data, Fourier transform-based data, and integral transform-based data (e.g., integral intensity, mean time of flight, time-of-flight variance, Laplace-Merlin transform data). Full-time-resolution data offers high parameter estimation accuracy but requires massive computation. For time-window data, early time windows exhibit lower deep spatial sensitivity, while later time windows offer greater deep spatial sensitivity; however, calculating the earlier time window information using the finite difference method still results in significant computational burden. Fourier transform-based data can effectively characterize photon flight distribution through multi-frequency signals, but the computational burden is still higher than some integral transform-based data due to the need to solve complex finite element equations. Integral transform data, on the other hand, maintains a certain level of deep spatial sensitivity while significantly reducing computational costs. However, existing techniques often use a single data type for reconstruction, either incurring high computational costs in pursuing theoretical accuracy while failing to improve spatial computational accuracy, or sacrificing deep spatial sensitivity in an attempt to simplify computation, thus limiting the overall performance of DOT algorithms. Summary of the Invention

[0006] To address the shortcomings of existing technologies, the purpose of this application is to provide a time-domain diffuse optical tomography reconstruction method and system, aiming to solve the problem that when time-domain diffuse optical tomography uses integral transform data for optical reconstruction, it often uses a single feature data type and sacrifices the sensitivity of deep space to simplify calculations, resulting in poor reliability and stability of diffuse optical tomography.

[0007] The first aspect of this application relates to a method for time-domain diffuse optical tomography reconstruction, specifically including the following steps: Step S1: Establish a three-dimensional geometric model of the biological tissue to be reconstructed, use the time-domain diffusion equation to describe the light transmission process in the biological tissue, and obtain the time-point diffusion function on the surface of the three-dimensional geometric model; extract feature data based on the time-point diffusion function, and draw the spatial sensitivity distribution map of the feature data respectively. Step S2: Using the mean square error between the feature data obtained in step S1 and the feature data of the time point diffusion function simulated based on the reconstruction results as the data fidelity term, and combining it with the regularization term to construct the objective function, thus constructing the inverse problem optimization model; Step S3: Based on the spatial sensitivity distribution map and node depth, select the feature data with high shallow spatial sensitivity and input them into the inverse problem optimization model to reconstruct the optical parameters of the shallow nodes; Step S4: Given that the optical parameters of the shallow nodes are known, select the feature data with high spatial sensitivity of the deep nodes based on the spatial sensitivity distribution map, and input them into the inverse problem optimization model to reconstruct the optical parameters of the deep nodes.

[0008] In some implementations, the time-domain diffusion equation is: ; in, Represents the absorption coefficient in biological tissues; Represents the diffusion coefficient; The reduced scattering coefficient; Represents the speed of light; Represents light intensity; Represents the light source item; This is the current position in the solution process; Time point spread function for: ; in, Represents the boundary normal derivative; It is the solution domain.

[0009] In some implementations, the characteristic data includes logarithmic integral intensity and mean flight time; The method for drawing a spatial sensitivity distribution map based on feature data in step S1 specifically includes the following steps: The time-point diffusion function, which varies with time, is collected at different times to obtain the logarithmic integral intensity and average flight time. The Jacobian matrix of the logarithmic integral intensity and the Jacobian matrix of the mean flight time are calculated based on the adjoint method. Combined with the sign function, the spatial sensitivity distribution maps of the logarithmic integral intensity and the mean flight time are obtained.

[0010] In some implementations, the zeroth and first moments of the time-point spread function are respectively: ; ; Logarithmic Integral Strength and average flight time for: ; ; in, Diffusion function at time points With respect to the zeroth moment of time; Diffusion function at time points The first moment with respect to time.

[0011] In some implementations, in step S3, a sensitivity threshold is set, and based on the spatial sensitivity distribution of the logarithmic integral intensity, nodes with spatial sensitivity higher than the sensitivity threshold are reconstructed using the logarithmic integral intensity, while nodes with spatial sensitivity lower than the sensitivity threshold are reconstructed using the average time of flight in step S4. Alternatively, a depth threshold for biological tissue can be set, and nodes whose depth from the detection surface is less than the depth threshold are reconstructed using logarithmic integral intensity; otherwise, optical parameters are reconstructed using the average time of flight in step S4. In step S3, when reconstructing optical parameters using logarithmic integral intensity, the absorption coefficient and diffusion coefficient are considered to be homogeneous. In step S4, when reconstructing optical parameters using the average time of flight, the absorption coefficient is considered heterogeneous and the diffusion coefficient is considered homogeneous.

[0012] In some implementations, the quasi-Newton method is used to solve the objective function in steps S3 and S4. Specifically, the quasi-Newton method is to use the adjoint method to solve the gradient of the objective function, and rewrite the product of the Jacobian matrix and the mean square error in the gradient expression of the objective function into the form of the product of the adjoint solution, the basis matrix, and the solution of the original problem and the mean square error through the associative law of matrix multiplication.

[0013] In some implementations, the iteration termination condition of the quasi-Newton method is: the change in optical parameters between two adjacent iterations is less than a preset threshold, or the decrease in the objective function value is less than a preset threshold, or the preset maximum number of iterations is reached.

[0014] The second aspect of this application relates to a time-domain diffuse optical tomography reconstruction system, comprising: The time-point diffusion function acquisition module is used to establish a three-dimensional geometric model of the biological tissue to be reconstructed, use the time-domain diffusion equation to describe the light transmission process in the biological tissue, and obtain the time-point diffusion function on the surface of the three-dimensional geometric model. The sensitivity distribution plotting module is used to extract feature data based on the time-point spread function and plot the spatial sensitivity distribution of the feature data. The optimization model building module is used to construct an objective function by using the mean square error between the feature data extracted from the sensitivity distribution plotting module and the feature data of the diffusion function at the time point simulated based on the reconstruction results as the data fidelity term, and using the regularization term as the output to establish an inverse problem optimization model. The shallow node reconstruction module is used to select feature data with high shallow spatial sensitivity based on the spatial sensitivity distribution map and input them into the inverse problem optimization model to reconstruct the optical parameters of the shallow nodes. The deep node reconstruction module is used to reconstruct the optical parameters of deep nodes by selecting feature data with high spatial sensitivity based on the spatial sensitivity distribution map, given that the optical parameters of shallow nodes are known, and inputting them into the inverse problem optimization model.

[0015] In some implementations, the feature data includes logarithmic integral intensity and mean flight time; the sensitivity distribution plotting module includes a feature extraction unit and a distribution plotting unit; The feature extraction unit is used to collect the time-point spread function at different times to obtain the logarithmic integral intensity and average flight time. The distribution map plotting unit is used to calculate the Jacobian matrix of the logarithmic integral intensity and the Jacobian matrix of the mean flight time based on the adjoint method, and combined with the sign function, to obtain the spatial sensitivity distribution map of the logarithmic integral intensity and the mean flight time.

[0016] In some implementations, the shallow node reconstruction module and the deep node reconstruction module use the quasi-Newton method to solve the objective function. Specifically, the quasi-Newton method is to use the adjoint method to solve the gradient of the objective function, and through the matrix multiplication associative law, rewrite the product of the Jacobian matrix and the mean square error in the gradient expression of the objective function into the form of the adjoint problem solution, the basis matrix, and the original problem solution and the mean square error.

[0017] In some implementations, the iteration termination condition of the quasi-Newton method in the shallow node reconstruction module and the deep node reconstruction module is: the change in optical parameters between two adjacent iterations is less than a preset threshold, or the decrease in the objective function value is less than a preset threshold, or the preset maximum number of iterations is reached.

[0018] In some implementations, the zeroth and first moments of the time-point spread function are respectively: ; ; Logarithmic Integral Strength and average flight time for: ; ; in, Diffusion function at time points With respect to the zeroth moment of time; Diffusion function at time points The first moment with respect to time.

[0019] In some implementations, the optical parameters include the absorption coefficient and the diffusion coefficient; the time-domain diffusion equation in the time-point diffusion function acquisition module is: ; in, Represents the absorption coefficient in biological tissues; Represents the diffusion coefficient; The reduced scattering coefficient; Represents the speed of light; Represents light intensity; Represents the light source item; This is the current position in the solution process; Time point spread function for: ; in, Represents the boundary normal derivative; It is the solution domain.

[0020] Thirdly, this application provides an electronic device, comprising: at least one memory for storing a program; and at least one processor for executing the program stored in the memory, wherein when the program stored in the memory is executed, the processor is configured to execute the method described in the first aspect or any possible implementation thereof.

[0021] Fourthly, this application provides a computer-readable storage medium storing a computer program that, when run on a processor, causes the processor to perform the method described in the first aspect or any possible implementation thereof.

[0022] Fifthly, this application provides a computer program product that, when run on a processor, causes the processor to perform the method described in the first aspect or any possible implementation thereof.

[0023] It is understood that the beneficial effects of the second to fifth aspects mentioned above can be found in the relevant descriptions in the first aspect mentioned above, and will not be repeated here.

[0024] Overall, the technical solutions conceived in this application have the following beneficial effects compared with the prior art: This application provides a temporal diffuse optical tomography reconstruction method. Based on spatial sensitivity distribution and depth information, it determines the detection depth of each feature data. More specifically, the logarithmic integral intensity is more sensitive to changes in the absorption coefficient of shallow tissues, while the mean time of flight is more sensitive to the absorption coefficient of deep tissues. This provides a basis for subsequent reconstruction of shallow and deep tissues based on a partitioning strategy, improving the stability and reliability of the reconstruction results. Furthermore, this method can provide reliable technical support for depth sensing in the field of optical imaging. This application uses logarithmic integral intensity and mean time of flight to reconstruct shallow and deep tissues respectively. During the reconstruction process, non-target regions are treated as homogeneous dimensionality reduction. Therefore, compared with existing methods that use single data features for reconstruction, the spatial sensitivity of the data corresponding to the reconstruction nodes is improved, and the space of the optical parameters to be reconstructed is reduced, resulting in a reconstruction accuracy of less than 15% for optical parameters.

[0025] This application provides a temporal diffuse optical tomography reconstruction method that improves the gradient calculation process in the quasi-Newton method. The product of the Jacobian matrix and the mean square error in the gradient expression of the objective function is rewritten as the product of the adjoint problem solution, the basis matrix, the original problem solution, and the mean square error. Compared with existing reconstruction methods based on Newton's method, it does not require explicit solution and storage of all Jacobian matrices, which significantly reduces memory usage and computational overhead, and effectively improves the algorithm's running efficiency by at least 50%.

[0026] This application provides a temporal-domain diffuse optical tomography reconstruction method that differentiates the processing for each tissue layer. Specifically, a sensitivity threshold is set, and based on the spatial sensitivity distribution of logarithmic integral intensity, nodes with spatial sensitivity higher than the threshold are reconstructed using logarithmic integral intensity, while nodes with spatial sensitivity lower than the threshold are reconstructed using average time-of-flight. Alternatively, a depth threshold is set for biological tissues; nodes with a depth less than the detection surface are reconstructed using logarithmic integral intensity, otherwise, average time-of-flight is used to reconstruct optical parameters. In step S3, when reconstructing the optical parameters of superficial tissues (skull and scalp) using logarithmic integral intensity, both absorption and diffusion coefficients are considered homogeneous. However, in step S4, the absorption coefficient is considered heterogeneous, while the diffusion coefficient is considered homogeneous. This processing significantly reduces the dimensionality of unknowns in the inverse problem, effectively alleviating the pathological problem. In step S4, the optical parameters reconstructed from superficial tissues are used as known quantities to eliminate the interference of surface optical parameters on the inversion process of deep parameters, maintaining the spatial calculation accuracy of the target tissue layer. Attached Figure Description

[0027] Figure 1 This is a schematic flowchart of the time-domain diffuse optical tomography reconstruction method provided in the embodiments of this application.

[0028] Figure 2 These are MRI cross-sectional images of the human brain model provided in the embodiments of this application in the coronal, sagittal and transverse planes.

[0029] Figure 3 This is a schematic diagram of the source-probe arrangement on the human brain model provided in the embodiments of this application.

[0030] Figure 4 This is a set of spatial sensitivity distribution maps for a specific source-detector pair provided in the embodiments of this application.

[0031] Figure 5 This is a spatial distribution map of the difference between the normalized average flight time and the normalized logarithmic integral intensity sensitivity provided in the embodiments of this application.

[0032] Figure 6 This is a layered reconstruction result diagram of the absorption coefficient provided in the embodiments of this application. Detailed Implementation

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

[0034] In this application, the term "and / or" describes the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent three cases: A existing alone, A and B existing simultaneously, and B existing alone. In this application, the symbol " / " indicates that the related objects are in an "or" relationship, for example, A / B means A or B.

[0035] In this application, the terms “first” and “second” are used to distinguish different objects, rather than to describe a specific order of objects.

[0036] In the embodiments of this application, the terms "exemplary" or "for example" are used to indicate that something is an example, illustration, or description. Any embodiment or design that is described as "exemplary" or "for example" in the embodiments of this application should not be construed as being more preferred or advantageous than other embodiments or design. Specifically, the use of the terms "exemplary" or "for example" is intended to present the relevant concepts in a specific manner.

[0037] In the description of the embodiments in this application, unless otherwise stated, "multiple" means two or more.

[0038] The embodiments of this application are described below with reference to the accompanying drawings.

[0039] like Figure 1 As shown, this application provides a time-domain diffuse optical tomography reconstruction method, including the following steps: Step S1: Establish a three-dimensional geometric model of the tissue to be reconstructed. Perform tetrahedral meshing on the three-dimensional geometric model to generate discrete nodes and elements. Label different tissue nodes according to the actual anatomical structure. Use the time-domain diffusion equation to describe the light transmission process in biological tissue. Assign different optical parameters to each node. Solve the time-domain diffusion equation using the finite element method to obtain the time-point diffusion function on the surface of the three-dimensional geometric model, which serves as the reference data for solving the subsequent inverse problem. The optical parameters include the absorption coefficient and diffusion coefficient in the diffusion equation. The scattering process of light in biological tissue is based on the time-domain diffusion equation. In order for the three-dimensional geometric model to conform to the optical characteristics of actual tissue, the absorption coefficient must be much smaller than the scattering coefficient.

[0040] Feature data is extracted based on the time-point spread function, and spatial sensitivity distribution maps of the feature data are plotted to serve as the basis for allocating feature data for reconstruction to each node; the feature data includes logarithmic integral intensity and mean flight time. Step S2: The mean square error between the feature data of the time point diffusion function obtained in step S1 and the feature data of the time point diffusion function obtained by forward simulation based on the reconstructed optical parameters is used as the data fidelity term. Combined with the regularization term, the objective function is constructed to build the inverse problem optimization model. Step S3: Select feature data with high sensitivity in shallow space, and reconstruct the optical parameters of shallow nodes according to the inverse problem optimization model in step S2. The optical parameters include absorption coefficient and diffusion coefficient. Assign these optical parameters to the shallower nodes in step S1. Step S4: Given that the optical parameters of the shallower tissue obtained from the inversion in step S3 are known, the optical parameters of the deeper nodes are reconstructed according to the inverse problem optimization model in step S2. The optical parameters include the absorption coefficient and the diffusion coefficient until the reconstruction is completed.

[0041] In steps S3 and S4, depending on whether the layer of tissue belongs to the target tissue, if the layer of tissue belongs to the target tissue, the optical parameters to be reconstructed are set to non-uniform; otherwise, they are set to uniform. This reduces the dimensionality of the inverse problem while controlling the discrete error of the target tissue layer.

[0042] In some implementations, the method for solving the objective function using the quasi-Newton method in step S2 is as follows: the gradient calculation involved in the Jacobian matrix of the feature data is calculated using the adjoint method, and the computational complexity and memory usage are reduced by using the matrix multiplication associative law and matrix reconstruction techniques. Furthermore, the matrix method includes: calculating the Jacobian matrix of the logarithmic integral intensity and the Jacobian matrix of the average flight time based on the adjoint method, avoiding the explicit construction and storage of the complete Jacobian matrix; on this basis, using the associative law of matrix multiplication, the product of the Jacobian matrix and the mean square error vector in the gradient expression is rewritten as the product of the adjoint problem solution, the basis matrix, the original problem solution and the error vector, thereby significantly reducing the computational overhead and storage requirements, and effectively improving the overall running efficiency of the algorithm; Example 1 This application uses the logarithmic integral intensity and mean time of flight derived from the time-point spread function as examples, employing a human brain model and reflective detection geometry to demonstrate its advantages in the field of brain functional imaging. Specific details of the three-dimensional geometric model used in this application can be found in [link to application]. Figure 2 and Figure 3 Based on this three-dimensional geometric model, a three-layer model including the scalp, skull, and brain parenchyma was constructed. Step S1: Establish a three-dimensional geometric model of the tissue to be reconstructed, perform tetrahedral meshing on the three-dimensional geometric model to generate discrete nodes and elements, label different tissue nodes according to the real anatomical structure, use the time-domain diffusion equation to describe the light transmission process in biological tissue, assign different optical parameters to each node, use the finite element method to solve the time-domain diffusion equation, and obtain the time-point diffusion function on the surface of the three-dimensional geometric model. More specifically, the time-domain diffusion equation is: (1); in, Represents the absorption coefficient in biological tissues; Represents the diffusion coefficient; The reduced scattering coefficient; Represents the speed of light; Represents light intensity; The light source term represents the light source term; in this application, given the extremely short pulse width of the pulsed light source used, the light source term can be approximated as the Dirac function in spacetime. , The function is the Dirac function; the pulse start time is taken as the time zero point. This is the location of the current solution to the requirement; Position of the light source; For boundary conditions, this application uses Robin boundary conditions: (2); in, It is the reflection coefficient within the diffusion transmission; Represents the boundary normal derivative; It is the solution domain; It is the boundary normal vector; The diffusion function for the detected time points is given by the following equation: (3); Combining the Robin boundary conditions, we can obtain in, .

[0043] Within the finite element framework, we have: (4); in, and It can be represented as: (5); (6); in, k Representing the k Unit; The total number of units with a constant diffusion coefficient; For the first k Diffusion coefficient at each node; For the first k Diffusion coefficient per unit cell; This represents the total number of units with a constant absorption coefficient. Representativeness and diffusion coefficient The basis matrix of the relevant stiffness matrix; Representative and absorption coefficient The basis matrix of the relevant stiffness matrix; and The elements in are as follows: (7); (8); in, represent The first in One element, represent The first in One element, and Represents the basis functions for finite element interpolation; It is a symmetric matrix.

[0044] After obtaining the finite element solution vector, the time-point spread function is passed through the boundary measurement operator. Map the solution vector to the probe vector Elements in matrix B Representing the i The position vectors of each detector; The interpolation basis function represents the unit interpolation basis function, which is only used when... There is a unit value above; The purpose of step S1 is to build a model and perform physical simulation. In order to solve the model, tetrahedral meshing is first required. This application uses the TetGen interface provided by PyVista and an adaptive clustering algorithm to optimize the mesh, so as to enhance the stability of subsequent finite element calculations.

[0045] Feature data were extracted based on the time-point spread function, and spatial sensitivity distribution maps of the feature data were plotted to serve as the basis for allocating feature data for each node to participate in the reconstruction. The feature data included logarithmic integral intensity and mean flight time. The method for obtaining the spatial sensitivity of the feature data is as follows: Step S1.1: Time-point spread functions that vary with time can be collected at different times. From these, the logarithmic integral intensity and mean flight time are derived as two characteristic data points for reconstruction. The operator is: (9); Let the zeroth and first moments of the time-point spread function be: (10); (11); Logarithmic Integral Strength and average flight time The definition is as follows: (12); (13); in, Diffusion function at time points With respect to the zeroth moment of time; Diffusion function at time points With respect to the first moment of time; Step S1.2: The Jacobian matrix of the time-point spread function with respect to optical parameters can be divided into blocks according to the detection data obtained under different light source illumination. For the Jacobian submatrix of the logarithmic integral strength Among the elements According to the adjoint method, we have: (14); in, For the first i A solution to an adjoint problem; The basis matrix of the stiffness matrix can represent or ; For the first j The solution to the original problem; For the first The light source emits light, the first... The integral intensity signal of the spread function at each time point detected by the detector; This represents the zeroth moment with respect to time; For the Jacobian submatrix of average flight time Among them, elements According to the differentiation rule, we have: (15); in, For the first The optical parameters of each unit represent or ; The zeroth moment of the diffusion function at time points; The first moment of the diffusion function at time points; and The following calculation can be performed using the adjoint method: (16); (17); in, The first moment with respect to time; Because diffuse optical computed tomography (DCT) is limited by spatial sensitivity, its ability to detect deep tissues is limited. Therefore, it is necessary to determine the effective reconstruction depth based on the spatial sensitivity matrix. Within the effective detection range, different measurement data have different sensitivities to different tissue depths. The purpose of step S2 is to determine the depth detection range of the logarithmic integral intensity and the average time of flight, and to utilize the difference in spatial sensitivity between the two to achieve layered reconstruction of tissues at different depths. Step S1.3: Calculate the spatial sensitivity function, defined as follows: (18); Among them, the sign function Both the logarithmic integral intensity and the mean flight time have negative unit values.

[0046] Step S2: Using the mean square error between the reconstructed feature data and the feature data extracted based on the time-point diffusion function as input, construct the objective function by combining the regularization term, solve the objective function using the quasi-Newton method, and establish an inverse problem optimization model with the reconstructed optical parameters as output; wherein, the gradient calculation in the quasi-Newton method is based on the adjoint method, and the computational complexity and memory usage are reduced by matrix multiplication associative law and matrix reconstruction technology. More specifically, the optical parameter reconstruction method under the finite element framework proposed in step S2 establishes the following optimization objective: (19); in, The emission data is based on simulation or experimental measurements. The emission data is calculated based on modeling. It is a parameter vector composed of absorption coefficient and diffusion coefficient; This is the regularization term; the method can be solved using the quasi-Newton method. For the gradient in the quasi-Newton method, this application proposes a fast solution method. For the error term... The gradients are: (20); The Jacobian matrices for the logarithmic integral intensity and the mean flight time can be calculated according to formulas (14) and (15), respectively. For the integral strength Jacobian matrix The submatrix can be written using the adjoint method as follows: (twenty one); in, (twenty two); (twenty three); (twenty four); Its gradient is: (25); Considering the high detector reuse rate and the fact that the number of detectors is much smaller than the total number of nodes after tetrahedral partitioning, and Combining calculations with previous terms reduces memory usage and saves computational resources, and eliminates the need to store the entire Jacobian matrix.

[0047] For the Jacobian matrix of average flight time The submatrix can be written using the adjoint method as follows: (26); in, Represents matrix dot product, and its gradient is: (27); Similarly, in order to conserve resources, priority should be given to... The zeroth and first moments and Combine calculations.

[0048] Step S3: Select feature data with high sensitivity in shallow space, and reconstruct the optical parameters of shallow nodes according to the inverse problem optimization model in Step S2. The optical parameters include absorption coefficient and diffusion coefficient. Assign the optical parameters to the shallower nodes in Step S2. Step S4: Given that the optical parameters of the shallower tissue obtained from the inversion in step S3 are known, the optical parameters of the deeper nodes are reconstructed according to the inverse problem optimization model in step S2. The optical parameters include the absorption coefficient and the diffusion coefficient until the reconstruction is completed.

[0049] In near-infrared brain functional imaging, the dynamic changes in hemoglobin concentration are related to the absorption coefficient of brain parenchyma. Therefore, this application considers the absorption and diffusion coefficients of the scalp and skull as homogeneous, and the diffusion coefficient of brain parenchyma as homogeneous, while only considering its absorption coefficient as heterogeneous. This processing method can effectively reduce the dimensionality of the parameter space while controlling the discrete error of brain parenchyma.

[0050] Based on the initial optical parameter estimates, the spatial sensitivity function is calculated according to step S2, and its spatial distribution is plotted. Figure 4 The spatial sensitivity distribution of the absorption coefficient to the mean flight time is shown at a specific source-detector distance. Figure 5This demonstrates the difference in spatial sensitivity distribution between the normalized logarithmic integral intensity and the mean time of flight to the absorption coefficient, intuitively reflecting their response characteristics at different depths. The results show that the logarithmic integral intensity is more sensitive to changes in the absorption coefficient of shallow tissues, while the mean time of flight is more sensitive to changes in the absorption coefficient of deep tissues. This provides a basis for subsequent reconstruction of shallow and deep tissues based on a partitioning strategy. In the simulation calculation, the sensitivity information corresponding to all source-detector pairs is integrated to classify discrete nodes into shallow and deep regions. In this numerical experiment, the scalp and skull are classified as shallow tissues, and the brain parenchyma is defined as deep tissue. In step S3, the absorption and diffusion coefficients of superficial tissues (i.e., skull and scalp) are reconstructed using the logarithmic integral intensity, and these are treated as homogeneous parameters in a homogeneous medium (both absorption and diffusion coefficients are considered homogeneous). This approach significantly reduces the dimensionality of unknowns in the inverse problem, thereby effectively alleviating the ill-conditioned problem. Using the mean square error between the calculated value of the logarithmic integral intensity and the feature data extracted from the diffusion function at time points as the objective function, based on the inverse problem solution framework established in step S3, iterative optimization is performed using a quasi-Newton method without introducing a regularization term, until a preset iteration termination condition is met, obtaining the optical parameters of the tissue. The above results are then substituted into the subsequent deep tissue reconstruction model as known quantities to eliminate the interference of surface optical parameters on the deep parameter inversion process.

[0051] Step S4: Given the surface optical parameters obtained in Step S3, a heterogeneous parameter inversion model for deep tissue is constructed using only the optical parameters of deep nodes as unknowns. The mean square error between the calculated average time of flight and the feature data extracted from the time-point diffusion function is used as the objective function. Based on the inverse problem-solving framework established in Step S2, a regularization term is introduced to suppress ill-posedness in deep reconstruction. Iterative optimization is performed using a quasi-Newton method until a preset iteration termination condition is met. During the iteration process, the surface parameters remain unchanged, and only the deep node parameters are updated. The absorption coefficient is considered heterogeneous, and the diffusion coefficient is considered homogeneous. Finally, the spatial distribution of the deep tissue optical parameters is obtained, completing the entire process of layered reconstruction. Figure 6 The reconstruction results are shown using the absorption coefficient as an example.

[0052] Example 2 The second aspect of this application relates to a time-domain diffuse optical tomography reconstruction system, comprising: The time-point diffusion function acquisition module is used to establish a three-dimensional geometric model of the biological tissue to be reconstructed, use the time-domain diffusion equation to describe the light transmission process in the biological tissue, and obtain the time-point diffusion function on the surface of the three-dimensional geometric model. The sensitivity distribution plotting module is used to extract feature data based on the time-point spread function and plot the spatial sensitivity distribution of the feature data. The optimization model building module is used to construct an objective function by using the mean square error between the feature data extracted from the sensitivity distribution plotting module and the feature data of the diffusion function at the time point simulated based on the reconstruction results as the data fidelity term, and using the regularization term as the output to establish an inverse problem optimization model. The shallow node reconstruction module is used to select feature data with high shallow spatial sensitivity based on the spatial sensitivity distribution map and input them into the inverse problem optimization model to reconstruct the optical parameters of the shallow nodes. The deep node reconstruction module is used to reconstruct the optical parameters of deep nodes by selecting feature data with high spatial sensitivity based on the spatial sensitivity distribution map, given that the optical parameters of shallow nodes are known, and inputting them into the inverse problem optimization model.

[0053] In some implementations, the feature data includes logarithmic integral intensity and mean flight time; the sensitivity distribution plotting module includes a feature extraction unit and a distribution plotting unit; The feature extraction unit is used to collect the time-point spread function at different times to obtain the logarithmic integral intensity and average flight time. The distribution map plotting unit is used to calculate the Jacobian matrix of the logarithmic integral intensity and the Jacobian matrix of the mean flight time based on the adjoint method, and combined with the sign function, to obtain the spatial sensitivity distribution map of the logarithmic integral intensity and the mean flight time.

[0054] In some implementations, the shallow node reconstruction module and the deep node reconstruction module use the quasi-Newton method to solve the objective function. Specifically, the quasi-Newton method is to use the adjoint method to solve the gradient of the objective function, and through the matrix multiplication associative law, rewrite the product of the Jacobian matrix and the mean square error in the gradient expression of the objective function into the form of the adjoint problem solution, the basis matrix, and the original problem solution and the mean square error.

[0055] In some implementations, the iteration termination condition of the quasi-Newton method in the shallow node reconstruction module and the deep node reconstruction module is: the change in optical parameters between two adjacent iterations is less than a preset threshold, or the decrease in the objective function value is less than a preset threshold, or the preset maximum number of iterations is reached.

[0056] In some implementations, the zeroth and first moments of the time-point spread function are respectively: ; ; Logarithmic Integral Strength and average flight time for: ; ; in, Diffusion function at time points With respect to the zeroth moment of time; Diffusion function at time points The first moment with respect to time.

[0057] In some implementations, the optical parameters include the absorption coefficient and the diffusion coefficient; the time-domain diffusion equation in the time-point diffusion function acquisition module is: ; in, Represents the absorption coefficient in biological tissues; Represents the diffusion coefficient; The reduced scattering coefficient; Represents the speed of light; Represents light intensity; Represents the light source item; This is the current position in the solution process; Time point spread function for: ; in, Represents the boundary normal derivative; It is the solution domain.

[0058] Compared with the prior art, this application has the following advantages: This application provides a temporal diffuse optical tomography reconstruction method. Based on spatial sensitivity distribution and depth information, it determines the detection depth of each feature data and accordingly determines the feature data used for reconstruction of each node. This method can provide reliable technical support for depth sensing in the field of optical imaging.

[0059] This application provides a temporal diffuse optical tomography reconstruction method, which significantly improves the gradient calculation process in the quasi-Newton method. It utilizes matrix operations to significantly reduce memory usage and computational overhead, effectively improving the overall algorithm's running efficiency.

[0060] This application provides a temporal diffuse optical tomography reconstruction method that employs a layered reconstruction approach, performing differentiated processing on each tissue layer. This effectively reduces the spatial dimensionality of parameters while maintaining the spatial computational accuracy of the target tissue layer. This method not only alleviates the ill-conditioned nature of the reconstruction problem but also effectively reduces computational overhead, significantly improving the stability and reliability of the reconstruction results.

[0061] It should be understood that the above-described device is used to execute the methods in the above embodiments. The implementation principle and technical effect of the corresponding program modules in the device are similar to those described in the above methods. The working process of the device can be referred to the corresponding process in the above methods, and will not be repeated here.

[0062] Based on the methods in the above embodiments, this application provides an electronic device that may include a processor, a communications interface, a memory, and a communication bus, wherein the processor, communications interface, and memory communicate with each other via the communication bus. The processor may invoke logical instructions stored in the memory to execute the methods in the above embodiments.

[0063] Furthermore, the logical instructions in the aforementioned memory can be implemented as software functional units and, when sold or used as independent products, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or a portion of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this application.

[0064] Based on the methods in the above embodiments, this application provides a computer-readable storage medium storing a computer program that, when run on a processor, causes the processor to execute the methods in the above embodiments.

[0065] Based on the methods in the above embodiments, this application provides a computer program product that, when run on a processor, causes the processor to execute the methods in the above embodiments.

[0066] It is understood that the processor in the embodiments of this application can be a central processing unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, transistor logic devices, hardware components, or any combination thereof. A general-purpose processor can be a microprocessor or any conventional processor.

[0067] The method steps in this application embodiment can be implemented in hardware or by a processor executing software instructions. The software instructions can consist of corresponding software modules, which can be stored in random access memory (RAM), flash memory, read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), electrically erasable programmable read-only memory (EEPROM), registers, hard disks, portable hard disks, CD-ROMs, or any other form of storage medium known in the art. An exemplary storage medium is coupled to the processor, enabling the processor to read information from and write information to the storage medium. Of course, the storage medium can also be a component of the processor. The processor and the storage medium can reside in an ASIC.

[0068] In the above embodiments, implementation can be achieved entirely or partially through software, hardware, firmware, or any combination thereof. When implemented using software, it can be implemented entirely or partially as a computer program product. The computer program product includes one or more computer instructions. When the computer program instructions are loaded and executed on a computer, all or part of the processes or functions described in the embodiments of this application are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted through the computer-readable storage medium. The computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (e.g., coaxial cable, fiber optic, digital subscriber line (DSL)) or wireless (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center that integrates one or more available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., DVD), or a semiconductor medium (e.g., solid-state disk (SSD)).

[0069] It is understood that the various numerical designations used in the embodiments of this application are merely for the convenience of description and are not intended to limit the scope of the embodiments of this application.

[0070] Those skilled in the art will readily understand that the above description is merely a preferred embodiment of this application and is not intended to limit this application. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this application should be included within the scope of protection of this application.

Claims

1. A method for reconstructing time-domain diffuse optical tomography, characterized in that, Includes the following steps: Step S1: Establish a three-dimensional geometric model of the biological tissue to be reconstructed, use the time-domain diffusion equation to describe the light transmission process in the biological tissue, and obtain the time-point diffusion function on the surface of the three-dimensional geometric model. Then, feature data are extracted based on the time-point spread function, and spatial sensitivity distribution maps of the feature data are plotted respectively. Step S2: Using the mean square error between the feature data obtained in step S1 and the feature data of the time point diffusion function simulated based on the reconstruction results as the data fidelity term, and combining it with the regularization term to construct the objective function, thus constructing the inverse problem optimization model; Step S3: Select feature data with high shallow spatial sensitivity based on the spatial sensitivity distribution map and input them into the inverse problem optimization model to reconstruct the optical parameters of the shallow nodes; Step S4: Given that the optical parameters of the shallow nodes are known, select the feature data with high spatial sensitivity of the deep nodes based on the spatial sensitivity distribution map, and input them into the inverse problem optimization model to reconstruct the optical parameters of the deep nodes. Among them, the characteristic data include logarithmic integral intensity and mean flight time; In step S3, a sensitivity threshold is set. Based on the spatial sensitivity distribution of the logarithmic integral intensity, nodes with spatial sensitivity higher than the sensitivity threshold are reconstructed using the logarithmic integral intensity. For nodes with spatial sensitivity lower than the sensitivity threshold, the optical parameters are reconstructed using the average time of flight in step S4. Alternatively, a depth threshold for biological tissue can be set, and nodes with a depth less than the detection surface can be reconstructed using logarithmic integral intensity reconstruction; otherwise, optical parameters can be reconstructed using average time of flight in step S4.

2. The temporal diffuse optical tomography reconstruction method according to claim 1, characterized in that, Optical parameters include absorption coefficient and diffusion coefficient; the time-domain diffusion equation is: ； in, Represents the absorption coefficient in biological tissues; Represents the diffusion coefficient; The reduced scattering coefficient; Represents the speed of light; Represents light intensity; Represents the light source item; This is the current position in the solution process; Time point spread function for: ; in, Represents the boundary normal derivative; It is the solution domain.

3. The time-domain diffuse optical tomography reconstruction method according to claim 1, characterized in that, The method for drawing a spatial sensitivity distribution map based on feature data in step S1 specifically includes the following steps: The time-point diffusion function, which varies with time, is collected at different times to obtain the logarithmic integral intensity and average flight time. The Jacobian matrix of the logarithmic integral intensity and the Jacobian matrix of the mean flight time are calculated based on the adjoint method, and the spatial sensitivity distribution maps of the logarithmic integral intensity and the mean flight time are calculated based on the Jacobian matrix.

4. The time-domain diffuse optical tomography reconstruction method according to claim 3, characterized in that, The zeroth and first moments of the time-point spread function with respect to time are as follows: ； ； Logarithmic Integral Strength and average flight time for: ； ； in, Diffusion function at time points With respect to the zeroth moment of time; Diffusion function at time points The first moment with respect to time.

5. The time-domain diffuse optical tomography reconstruction method according to claim 3 or 4, characterized in that, In step S3, when reconstructing optical parameters using logarithmic integral intensity, the absorption coefficient and diffusion coefficient are considered to be homogeneous. In step S4, when reconstructing optical parameters using the average time of flight, the absorption coefficient is considered heterogeneous and the diffusion coefficient is considered homogeneous.

6. The time-domain diffuse optical tomography reconstruction method according to claim 1, characterized in that, In steps S3 and S4, the quasi-Newton method is used to solve the objective function. Specifically, the quasi-Newton method is to use the adjoint method to solve the gradient of the objective function, and through the matrix multiplication associative law, rewrite the product of the Jacobian matrix and the mean square error in the gradient expression of the objective function into the form of the product of the adjoint problem solution, the basis matrix, and the original problem solution and the mean square error.

7. A time-domain diffuse optical tomography reconstruction system, characterized in that, include: The time-point diffusion function acquisition module is used to establish a three-dimensional geometric model of the biological tissue to be reconstructed, use the time-domain diffusion equation to describe the light transmission process in the biological tissue, and obtain the time-point diffusion function on the surface of the three-dimensional geometric model. The sensitivity distribution plotting module is used to extract feature data based on the time-point spread function and plot the spatial sensitivity distribution of the feature data. The optimization model building module is used to construct an objective function by using the mean square error between the feature data extracted from the sensitivity distribution plotting module and the feature data of the diffusion function at the time point simulated based on the reconstruction results as the data fidelity term, and using the regularization term as the output to establish an inverse problem optimization model. The shallow node reconstruction module is used to select feature data with high shallow spatial sensitivity based on the spatial sensitivity distribution map and input them into the inverse problem optimization model to reconstruct the optical parameters of the shallow nodes. The deep node reconstruction module is used to reconstruct the optical parameters of deep nodes by selecting feature data with high spatial sensitivity based on the spatial sensitivity distribution map, given that the optical parameters of shallow nodes are known, and inputting them into the inverse problem optimization model. Among them, the characteristic data include logarithmic integral intensity and mean flight time; A sensitivity threshold is set, and based on the spatial sensitivity distribution of the logarithmic integral intensity, nodes with spatial sensitivity higher than the sensitivity threshold are reconstructed using the logarithmic integral intensity, while nodes with spatial sensitivity lower than the sensitivity threshold are reconstructed using the average time of flight. Alternatively, a depth threshold for biological tissue can be set. Nodes with a depth less than the detection surface are reconstructed using logarithmic integral intensity; otherwise, optical parameters are reconstructed using the average time-of-flight.

8. The time-domain diffuse optical tomography reconstruction system according to claim 7, characterized in that, The sensitivity distribution map plotting module includes a feature extraction unit and a distribution map plotting unit; The feature extraction unit is used to collect the time-point spread function at different times to obtain the logarithmic integral intensity and average flight time. The distribution map plotting unit is used to calculate the Jacobian matrix of the logarithmic integral intensity and the Jacobian matrix of the mean flight time based on the adjoint method, and to calculate the spatial sensitivity distribution map of the logarithmic integral intensity and the mean flight time based on the Jacobian matrix.

9. The time-domain diffuse optical tomography reconstruction system according to claim 7, characterized in that, The shallow node reconstruction module and the deep node reconstruction module use the quasi-Newton method to solve the objective function. Specifically, the quasi-Newton method is to use the adjoint method to solve the gradient of the objective function, and through the matrix multiplication associative law, rewrite the product of the Jacobian matrix and the mean square error in the gradient expression of the objective function into the form of the product of the adjoint problem solution, the basis matrix, and the original problem solution and the mean square error.

Images