A multi-frame photoacoustic image reconstruction method based on optical flow alignment and depth feature fusion

Abstract

The application discloses a kind of multi-frame photoacoustic image reconstruction methods based on optical flow alignment and depth feature fusion, the method comprises the following steps: step S1, photoacoustic signal data acquisition;Step S2, photoacoustic tomography image reconstruction;Step S3, optical flow calculation and image alignment;Step S4, depth feature fusion network training.The application is through optical flow motion correction and potential space learning mechanism, dynamically integrates the spatio-temporal information and complementary features in the multi-frame image after alignment, adaptively fuses the complementary information in the multi-frame image and suppresses noise, finally realizes the reconstruction of high signal-to-noise ratio and high spatial resolution image to biological tissue.Experimental results show that the method can significantly improve the image quality, restore the image distortion and detail loss caused by motion and noise, and provide a new effective scheme for promoting the robust clinical application of photoacoustic imaging technology.

Description

Technical Field

[0001] This invention belongs to the fields of medical imaging and biomedical engineering, and relates to a reconstruction method based on the fusion of image information from the same tomographic section at different time points. Specifically, it relates to a multi-frame photoacoustic image reconstruction method based on optical flow alignment and depth feature fusion. Background Technology

[0002] Photoacoustic tomography (PAT), as an emerging non-invasive medical imaging technique, has attracted widespread attention in recent years and shown promising clinical application prospects. This technique combines the high contrast of optical imaging with the deep penetration capability of ultrasound imaging, providing a highly promising imaging method through the photoacoustic effect. Specifically, pulsed lasers excite biological tissue to generate ultrasound signals, which are then detected by an ultrasound transducer. The intensity and profile of the photoacoustic signal are closely related to the tissue's optical absorption, which reflects the tissue's physiological and pathological state. By reconstructing the photoacoustic signal, the distribution of radiation absorption within the tissue can be revealed, thus providing important information for medical diagnosis.

[0003] However, the advantages of this technology largely depend on high-quality signal acquisition and image reconstruction. In actual in vivo imaging, motion artifacts caused by physiological movements (such as breathing and heartbeat) become a major obstacle affecting image clarity and quantitative accuracy. These artifacts lead to blurred fine structures and distorted edges, which not only reduces the diagnostic value of the image but also poses a serious challenge to traditional reconstruction methods. Specifically, it raises four interconnected core technical issues: First, the dynamic mismatch problem, namely, the non-rigid deformation between frames caused by physiological motion, which makes it impossible to accurately align multiple frames in space; second, the signal-to-noise ratio bottleneck problem, namely, in a living environment, single-frame images are often contaminated by noise, especially in deep tissues or low signal-to-noise ratio scenes, making it difficult to effectively extract weak signals; third, the adaptive fusion and reconstruction problem, namely, how to adaptively fuse complementary information from aligned but still noisy and partially informative multi-frame images to recover or even surpass the anatomical details of a single frame; and fourth, the high-quality training data dependency problem, namely, existing deep learning-based methods usually rely on "low-quality-high-quality" image pairs for supervised training, while obtaining such perfectly paired image pairs in real living imaging is extremely difficult, limiting the practical application of the algorithm.

[0004] While methods such as Universal Back Projection (UBP), iterative reconstruction algorithms, and time-reversal methods are relatively mature in static scenes, they perform poorly when dealing with dynamic interferences and struggle to simultaneously address the synergistic challenges of motion compensation, noise suppression, and information fusion. To address these issues, deep learning-based image reconstruction methods have received widespread attention in recent years and have shown significant potential in improving image quality. Through multi-level feature learning and complex nonlinear mapping, these methods can effectively improve the accuracy and robustness of image reconstruction. For example, the extractor-attention-predictor network designed by Wang et al. significantly improved the quantitative accuracy of photoacoustic imaging; the three-path feature transform network proposed by Ma et al. optimized the reconstruction quality of the ring array detection system; and the Y-Net developed by Lan et al. achieved good results in improving the robustness of in vivo imaging. However, despite their excellent performance in single-frame image reconstruction, most of these methods still employ a "single-set input - single-image output" reconstruction paradigm, failing to fully utilize the redundant information of the same anatomical section at consecutive time points in live imaging. Therefore, when faced with significant inter-frame motion and noise, single-frame-based reconstruction methods still struggle to systematically address the four problems mentioned above, particularly in achieving adaptive fusion of multi-frame information, and they still rely on difficult-to-obtain paired data during the training phase. Summary of the Invention

[0005] To overcome the limitations of existing photoacoustic imaging methods in dynamic live scenes, such as motion mismatch, noise interference, single-frame reconstruction, and dependence on high-quality data, this invention provides a multi-frame photoacoustic image reconstruction method based on optical flow alignment and deep feature fusion. The method framework comprises two stages: The first stage employs optical flow for pixel-level motion estimation and frame alignment, thereby addressing the dynamic mismatch problem, compensating for non-rigid deformation caused by physiological motion, and providing a spatially consistent sequence for subsequent fusion. The second stage constructs a deep feature fusion network based on self-supervised learning, utilizing the spatiotemporal consistency of the aligned multi-frame sequence to address the signal-to-noise ratio bottleneck and the adaptive fusion reconstruction problem. The core of this network is an encoder-fusion-decoder architecture, where a pyramid encoder extracts multi-scale features from each frame layer by layer, and a variational autoencoder (VAE) models and regularizes the latent space, enhancing the network's generalization ability and robustness. In the fusion stage, a densely connected network (DenseNet) module fuses deep features from multiple frames, and finally, a U-Net-style decoder progressively upsamples and refines the fused features to reconstruct a high-quality image. During optimization, mean squared error (MSE) and structural similarity (SSIM) are combined to ensure pixel fidelity and structural integrity, and the VAE latent space is constrained by KL divergence. An adaptive weighting mechanism is also employed to balance the loss, ensuring stable and efficient training. This method can fuse features in the latent space, achieving noise suppression and information enhancement, ultimately reconstructing images with high signal-to-noise ratio and high spatial resolution. This significantly improves reconstruction quality, reduces reliance on high-quality labeled data, and enhances the applicability of the method in real-world live-body imaging scenarios. Through this two-stage processing framework, this invention fundamentally overcomes the limitations of traditional methods in dynamic imaging and effectively solves the four key problems mentioned above. Experimental results show that this method can significantly improve image quality, restore image distortion and detail loss caused by motion and noise, and provide a new and effective solution for promoting the robust clinical application of photoacoustic imaging technology.

[0006] The objective of this invention is achieved through the following technical solution:

[0007] A multi-frame photoacoustic image reconstruction method based on optical flow alignment and depth feature fusion includes the following steps:

[0008] Step S1: Acquisition of photoacoustic signal data:

[0009] Step S1-1: Obtain simulation signal: Use the k-wave toolbox to simulate and generate photoacoustic signal;

[0010] Step S2, Photoacoustic Tomography Image Reconstruction:

[0011] Step S2-1: Set imaging parameters: Set the following parameters during the imaging process: radius of the annular ultrasonic detector, number of ultrasonic transducers, detection angle, size of the imaging physical area, pixel grid size, sampling frequency of the data acquisition unit, sound wave propagation speed, and aperture size.

[0012] Step S2-2: Image reconstruction using DAS imaging method: The image sequence is reconstructed using the DAS imaging method. In the DAS imaging method, the relationship between the beamforming signal and the received signal is defined as follows:

[0013]

[0014] in, This refers to the number of ultrasonic transducers. Indicates position Beamforming signal at the location, Indicates the first Each transducer array element at time... Received signal, Represents from position The generated sound waves propagate to the first The time required for each transducer array element;

[0015] Step S3, Optical Flow Calculation and Image Alignment:

[0016] Step S3-1: Calculate the optical flow field: Select one frame from the initial image sequence as a reference image, and calculate the pixel-level motion field (optical flow field) between the remaining target images and the reference image based on the optical flow method. The specific steps are as follows:

[0017] Input reference image and target image We model the neighborhood of each pixel using a quadratic polynomial: by performing a polynomial expansion, we locally fit a quadratic polynomial for each pixel in both images, in the following form:

[0018]

[0019] in, The coordinates are in local coordinates. It is a symmetric matrix. It is a vector. It is a scalar. Indicates the transpose operation;

[0020] A coarse-to-fine estimation is performed on the constructed image pyramid, starting from the coarsest level. Assuming an initial displacement field, at each level, based on the principle that polynomial coefficients should remain invariant under motion, the displacement vector of each pixel is solved by solving a weighted least-squares system constructed from all pixels in the neighborhood. This process satisfies the constant brightness constraint:

[0021]

[0022] in, Reference image in global space coordinates and time pixel brightness, This represents the pixel brightness of the target image at the corresponding moved position. These are the global spatial coordinates of the pixels in the image. The time of the current image frame. These are the horizontal and vertical displacements of pixels between two frames, respectively.

[0023] The displacement field is continuously updated through iterative optimization to minimize the intensity difference between the two images. The objective function is:

[0024]

[0025] in, Reference image in pixels The brightness value; Target image position after offset The brightness value; The index representing the pixel; Find the displacement that minimizes the error ;

[0026] The estimated optical flow at this level is upsampled to a preset number of levels (the specific number of levels is fine-tuned based on the image alignment effect) as the initial value. The above optimization process is repeated until the original resolution is reached. Finally, a dense optical flow field is output. ;

[0027] Step S3-2, Image Alignment: The calculated optical flow field is applied to each target image to perform geometric transformation, thereby obtaining an image sequence that is strictly aligned with the reference image in space. The specific steps are as follows:

[0028] Obtaining the optical flow field Then, by analyzing the target image Perform geometric deformation to make it match the reference image. Spatial alignment, for the aligned image Each coordinate position in Find its corresponding sub-pixel coordinates in the source image based on the optical flow field. Using bilinear interpolation from The pixel intensity value at that location is calculated and assigned to... After traversing all pixels, the final result is the same as the reference image. Precisely registered aligned images ;

[0029] Step S4: Training the deep feature fusion network:

[0030] Step S4-1, Data Preparation: Merge the channels of the aligned target image and the reference image, and perform normalization processing to form the network input;

[0031] Step S4-2, Pyramid Encoder Training: Input the processed image into the pyramid encoder for multi-scale feature extraction to obtain a high-quality feature map. Optimize the encoder parameters through gradient backpropagation. The specific steps are as follows:

[0032] Image stack The input pyramid encoder extracts multi-scale image features using a convolutional pyramid structure. The encoder employs... Layered convolutional structures are used to extract multi-scale feature representations of images:

[0033]

[0034] in, Representing the Multi-scale feature maps output after convolutional layers. The input image is used as the feature map of layer 0 of the network. No. The layer contains convolutional and batch-normalized convolutional operations. The LeakyReLU activation function is used. The index of the convolutional layers used in the network is represented as follows:

[0035]

[0036] in, It is the first network The output feature map generated by the layer To output the number of feature channels, This represents the spatial size of the output feature map. Represents the high-dimensional, multi-scale image features ultimately extracted by the entire encoder;

[0037] Step S4-3, VAE Encoder Training: Input the feature map obtained in step S4-2 into the encoder of the variational autoencoder (VAE). Use a Gaussian distribution to model the posterior distribution of the latent space, learn the mean vector and standard deviation vector of each set of features, and use reparameterization techniques to achieve differentiability of the sampling process, which facilitates gradient optimization. The specific steps are as follows:

[0038] Gaussian distribution approximation: The VAE encoder learns the mean and log-variance parameters through a neural network.

[0039]

[0040]

[0041] in, and These are potential embedded variables. The mean and variance, Indicates that the encoder processes the input. The hidden representation obtained afterwards and For different fully connected layers; the posterior distribution defined by the encoder is:

[0042]

[0043] in, It is a latent embedded variable The conditional distribution, given the input The subsequent distribution;

[0044] Reparameterization technique: To make the sampling process differentiable, reparameterization techniques are used to generate latent embedding variables. :

[0045]

[0046] in, , This represents element-wise multiplication;

[0047] Step S4-4, Fusion Layer Network Training: The fusion layer network receives all latent embedding variables from the VAE encoder. Information is integrated through convolutional layers with dense connections and reflection-filled structures to output a unified, information-rich fusion latent representation. The specific steps are as follows:

[0048] Multiple latent variable concatenation: The input latent variables are concatenated along the channel dimension:

[0049]

[0050] in, It involves multiple potential embedded variables. The result after stitching along the channel dimension. This represents the total number of channels;

[0051] Densely connected converged network: adopts Feature fusion is performed using a densely connected layer structure.

[0052]

[0053]

[0054] in, This indicates the total number of network layers, with a value of 5. Represents the index of the current layer. For the input feature map, For the first The output feature map of the layer, and finally, the fused latent variable features are:

[0055]

[0056] in, It is the final fusion feature;

[0057] Latent variable features are compressed into a single latent variable: fused features are compressed into the target dimension through multi-layer convolution.

[0058]

[0059]

[0060] in, To compress the total number of layers in the network, Represents the index of the current compression layer. No. Layers are used for convolutional operations to compress features. After the first The intermediate features obtained after layer compression ultimately yield a single latent variable after fusion:

[0061]

[0062] in, The final output, a compressed single latent variable, is the output of the last layer. ;

[0063] Steps S4-5: U-Net Decoder Training: The decoder receives the fused latent variables and gradually reconstructs the spatial structure of the features through convolution, upsampling, and skip connections. It then utilizes the ReLU activation function to enhance non-linear expressive power, ultimately outputting a high-quality fused image. The specific steps are as follows:

[0064] The fused latent variables are reconstructed into the original image. The U-Net architecture is used to implement the mapping from the latent space to the image space. The specific process is as follows:

[0065]

[0066] in, For the U-Net decoder, upsampling is performed step by step through a symmetrical encoder-decoder path and skip connections, and the final output is a reconstructed image with pixel values ​​in the range of [0, 1].

[0067] Step S4-6, Loss function calculation:

[0068] loss function It consists of three parts: structural similarity loss, mean squared error loss, and KL divergence loss, where:

[0069] Structural similarity loss Defined as:

[0070]

[0071] in, Represents the original input image. This represents the image reconstructed by the decoder. Represents the structural similarity index;

[0072] Mean square error loss Defined as:

[0073]

[0074] in, Denotes the Frobenius norm. , , , These represent the number of images to be fused, the number of channels in the images, and the height and width, respectively.

[0075] KL divergence loss Defined as:

[0076]

[0077] in, and They represent the latent variables respectively. The mean and standard deviation of each dimension The dimension representing the potential space;

[0078] loss function Defined as:

[0079]

[0080]

[0081] in, For balance coefficient, For the total reconstruction loss, As weight;

[0082] Step S4-7, Optimization process: The Adam optimizer is used to jointly optimize the parameters of the pyramid encoder, VAE encoder, fusion layer network and U-Net decoder by minimizing the loss function.

[0083] Compared to existing technologies, this invention significantly improves the reconstruction quality, robustness, and practicality of dynamic photoacoustic imaging by constructing a multi-frame photoacoustic image reconstruction method based on optical flow alignment and depth feature fusion. Specifically, it achieves the following:

[0084] 1. Dynamic Motion Artifact Correction: This method achieves pixel-level estimation and compensation of physiological motion in living tissue through optical flow calculation, thereby generating spatially aligned image sequences. This approach fundamentally solves the inter-frame mismatch problem in dynamic imaging, providing a spatiotemporally consistent data foundation for subsequent high-quality fusion.

[0085] 2. Deep complementary fusion of multi-frame information: A self-supervised deep feature fusion network is introduced to adaptively learn and fuse the spatiotemporal features and complementary information of aligned multi-frame images in the latent space, breaking through the information limitations of single-frame reconstruction and synergistically enhancing the signal-to-noise ratio and spatial resolution of the image.

[0086] 3. Self-supervised learning enhances generalization ability: The entire network adopts a self-supervised training paradigm, eliminating the need for difficult-to-obtain paired training data and optimizing solely using multiple frames of collected images. This effectively reduces the method's dependence on specific datasets, enhancing its generalization and deployability in practical applications.

[0087] 4. Joint Optimization Ensures Overall Performance: This invention places the motion correction module and the deep feature fusion module in two separate stages, and optimizes each stage using a multi-objective loss function. In the motion correction stage, the optical flow alignment module ensures accurate image alignment, providing high-quality input for subsequent feature fusion. In the deep feature fusion stage, self-supervised learning is used to fuse features from multiple frames, optimizing the signal-to-noise ratio and spatial resolution of the images. This two-stage optimization strategy ensures synergy between the stages, thereby achieving overall optimal performance in terms of high fidelity, structural similarity, and feature separability in image reconstruction. Attached Figure Description

[0088] Figure 1 This is a schematic diagram of the photoacoustic signal acquisition process;

[0089] Figure 2 This is a schematic diagram of a beamforming signal;

[0090] Figure 3This is a schematic diagram of forward optical flow and image alignment. (a) is the forward optical flow calculation module, which is used to estimate the optical flow field between images. (b) is the image alignment and registration module, which spatially aligns the images based on the estimated optical flow.

[0091] Figure 4 The diagram shows the fusion network architecture. (a) is the pyramid feature encoding module, (b) is the latent variable encoding module composed of variational autoencoder (VAE), (c) is the latent variable fusion module based on DenseNet, and (d) is the module that decodes latent variables into the image domain using the U-Net structure.

[0092] Figure 5 This is a schematic diagram illustrating the loss function calculation.

[0093] Figure 6 The following are simulation results under normal full sampling conditions (512 ultrasonic transducers, 360-degree detection angle): (a) is the tomographic image at the first moment, (b) is the tomographic image at the second moment, (c) is the tomographic image at the third moment, and (d) is the fused image obtained by the method of this invention. Detailed Implementation

[0094] The technical solution of the present invention will be further described below with reference to the accompanying drawings, but it is not limited thereto. Any modifications or equivalent substitutions to the technical solution of the present invention that do not depart from the spirit and scope of the technical solution of the present invention should be covered within the protection scope of the present invention.

[0095] This invention provides a multi-frame photoacoustic image reconstruction method based on optical flow alignment and depth feature fusion. The method dynamically integrates spatiotemporal information and complementary features from aligned multi-frame images through optical flow motion correction and latent spatial learning mechanisms. It adaptively fuses complementary information from multiple frames and suppresses noise, ultimately achieving high signal-to-noise ratio and high spatial resolution image reconstruction of biological tissues. Figures 1-5 As shown, the overall process is as follows:

[0096] Step S1: Acquisition of photoacoustic signal data:

[0097] Step S1-1: Obtain simulation signal: Use the k-wave toolbox to simulate and generate photoacoustic signal.

[0098] Step S2, Photoacoustic Tomography Image Reconstruction:

[0099] Step S2-1: Set imaging parameters: Set the parameters for the annular ultrasonic detector radius, number of ultrasonic transducers, detection angle, size of the imaging physical area, pixel grid size, data acquisition sampling frequency, sound wave propagation speed, and aperture size during the imaging process.

[0100] Step S2-2: Image reconstruction using DAS imaging method: Delay-and-Sum (DAS) imaging technology is applied to reconstruct the image sequence.

[0101] Step S3, Optical Flow Calculation and Image Alignment:

[0102] Step S3-1: Calculate the optical flow field: Select a frame from the initial image sequence as a reference image, and calculate the pixel-level motion field (optical flow field) between the remaining target images and the reference image based on the optical flow method.

[0103] Step S3-2, Image Alignment: Apply the calculated optical flow field to each target image and perform geometric transformation to obtain an image sequence that is strictly aligned with the reference image in space.

[0104] Step S4: Training the deep feature fusion network:

[0105] In this stage, adaptive feature learning and fusion are performed on the aligned image sequence through a self-supervised deep feature fusion network.

[0106] Step S4-1, Data Preparation: Merge the channels of the aligned target image and the reference image, and perform normalization processing to form the network input.

[0107] Step S4-2, Pyramid Encoder Training: The processed image is input into the pyramid encoder for multi-scale feature extraction to obtain a high-quality feature map. The encoder parameters are optimized through gradient backpropagation to improve its feature representation capability.

[0108] Step S4-3, VAE Encoder Training: Input the above feature maps into the encoder of the Variational Autoencoder (VAE), model the posterior distribution of the latent space using a Gaussian distribution, and learn the mean vector of each set of features. with standard deviation vector The reparameterization technique is used to achieve differentiability of the sampling process, which facilitates gradient optimization.

[0109] Step S4-4, Fusion Layer Network Training: The fusion layer network receives all latent embedding variables from the VAE encoder. Information is integrated through convolutional layers with dense connections and reflection-filled structures, resulting in a unified and information-rich fusion latent representation.

[0110] Steps S4-5: U-Net Decoder Training: The decoder receives the fused latent variables, gradually reconstructs the spatial structure of the features through convolution, upsampling, and skip connections, enhances the nonlinear expressive power using the ReLU activation function, and finally outputs a high-quality fused image.

[0111] Steps S4-6: Loss Function Calculation: The loss function consists of three parts: 1. Weighted structural similarity loss, used to preserve image structural details; 2. Weighted mean squared error loss, used to maintain intensity consistency; 3. KL divergence loss, used to constrain the latent spatial distribution to approximate a standard normal distribution, ensuring the separability and consistency of its features. The loss weights are adaptively determined by the amount of feature information extracted by the pre-trained VGG-16 network.

[0112] Steps S4-7: Optimization Process: The Adam optimizer is used to jointly optimize the parameters of the pyramid encoder, VAE encoder, fusion layer network, and U-Net decoder by minimizing the loss function. The training process includes forward propagation, loss calculation, and gradient backpropagation.

[0113] Compared with the prior art, the present invention has the following advantages:

[0114] 1. Effective motion artifact correction: By aligning with the image through optical flow calculation, the physiological motion of living tissue is explicitly estimated and compensated, fundamentally alleviating the motion mismatch problem in dynamic imaging and providing a solid foundation for high-quality image fusion.

[0115] 2. Excellent noise suppression capability: Adaptive feature learning and fusion in the latent space can effectively separate real tissue signals from imaging noise, significantly improve the signal-to-noise ratio of the output image, and enhance the robustness of the method in low-dose or high-noise scenarios.

[0116] 3. Overcoming the limitations of single-frame information: Innovatively utilizing the spatiotemporal complementary information between multi-frame aligned images, and adaptively fusing them through a deep learning network, the spatial resolution and detail fidelity of the images are comprehensively improved, overcoming the information bottleneck in traditional single-frame reconstruction.

[0117] 4. Reduced dependence on paired data: By adopting a self-supervised learning paradigm, network training does not require paired "high-quality-low-quality" image data. Optimization can be performed using only the collected multiple frames of images, improving the practicality and ease of use of the method.

[0118] 5. Self-supervised differentiable optimization: The entire process is constructed as a complete differentiable computation graph, supporting gradient propagation and parameter optimization from input to output. This ensures the collaborative work of the motion correction module and the feature fusion module under the self-supervised framework, achieving optimal overall performance.

[0119] Example:

[0120] Step 1: Signal Acquisition.

[0121] Simulation Signal Creation: The forward simulation of this photoacoustic imaging was conducted within a two-dimensional computational domain of approximately 11.7cm × 11.7cm, discretized into 900 × 900 grid points with a grid resolution of 130μm. The simulation used an initial sound pressure distribution extracted and normalized from a blood vessel morphology image as the sound source, embedded in a sound velocity field. The baseline sound velocity was set to 1512 m / s, and variations in the non-uniform distribution could be added based on external data. The density of the background medium was uniformly set to 1000 kg / m³. Sound wave propagation was numerically solved using the k-space pseudospectral method, with the time step automatically determined based on the sound velocity stability condition. The entire process involved 6000 time steps, and 20 perfectly matched layers were placed at the outer edge of the computational domain to absorb boundary reflections. Ultimately, the propagated sound waves were received by a complete ring-shaped ultrasonic sensor array containing 512 elements and a radius of 50 mm. The radio frequency signals of each element were recorded over time at a high sampling frequency of 80 MHz, thus completing the simulation of the complete physical process from photoacoustic excitation to ultrasonic signal acquisition.

[0122] Step 2: Signal to image reconstruction.

[0123] Imaging parameters: Circular ultrasonic transducer radius is 50 × 10⁻⁶. -3 The instrument has 512 ultrasonic transducers, a detection angle range of 360 degrees, an imaging area of ​​66.56mm × 66.56mm, 512 × 512 pixels, a sampling frequency of 80MHz, a sound velocity of 1496m / s, and an aperture size of 0.8mm.

[0124] Delayed Summation (DAS) Image Reconstruction: In the DAS algorithm, the relationship between the beamforming signal and the received signal can be defined as:

[0125]

[0126] in, This refers to the number of ultrasonic transducers. Indicates position Beamforming signal at the location, Indicates the first Each transducer array element at time... Received signal. Delay time. Represents from position The generated sound waves propagate to the first The time required for each transducer element. This delay time can be expressed as:

[0127]

[0128] in, Point With the The distance between each transducer array element The speed of sound.

[0129] Step 3: Optical flow calculation and image alignment.

[0130] Optical flow alignment: Input reference image and target image (All images are grayscale). First, a quadratic polynomial is used to model the neighborhood of each pixel. Specifically, a polynomial expansion is first performed to locally fit a quadratic polynomial for each pixel in both images, in the following form:

[0131]

[0132] in, Let A be the coordinates in local coordinates, A be a symmetric matrix used to capture local second-order changes, and b be a vector representing local first-order changes (gradients). It is a scalar, representing the bias term of the local signal (corresponding to the average value near the neighborhood center). This indicates the transpose operation. These coefficients are estimated using a weighted least-squares fit of the signal values ​​within the neighborhood. Subsequently, coarse-to-fine estimations are performed on the constructed image pyramid. Starting from the coarsest level, an initial displacement field (typically zero) is assumed. At each level, based on the principle that polynomial coefficients should remain invariant under motion, the displacement vector of each pixel is solved by solving a weighted least-squares system constructed from all pixels in the neighborhood. The core of this process is satisfying the constant brightness constraint:

[0133]

[0134] in, Reference image in global space coordinates and time pixel brightness, This represents the pixel brightness of the target image at the corresponding moved position. These are the global spatial coordinates of the pixels in the image. The time of the current image frame. These represent the horizontal and vertical displacements of pixels between two frames, respectively. The displacement field is continuously updated through iterative optimization to minimize the intensity difference between the two images. The objective function is:

[0135]

[0136] in, Reference image in pixels The brightness value; Target image position after offset The brightness value; The index representing the pixel; Find the displacement that minimizes the error Finally, the optical flow estimated at this level is upsampled to a finer level as the initial value, and the above optimization process is repeated until the original resolution is reached. Ultimately, a dense optical flow field is output. .

[0137] Image alignment: Obtaining optical flow field Then, by analyzing the target image Perform geometric deformation to make it match the reference image. Spatial alignment is achieved using inverse mapping and bilinear interpolation: for the aligned image Each coordinate position in Find its corresponding sub-pixel coordinates in the source image based on the optical flow field. .because Typically non-integer, therefore bilinear interpolation is used from The pixel intensity value at that location is calculated and assigned to... After traversing all pixels, the final result is a copy of the reference image. Precisely registered aligned images .

[0138] Step 4: Data preprocessing.

[0139] After normalizing the reference image and the aligned image, they are merged to form an image stack. .

[0140] Step 5: Self-supervised deep feature fusion.

[0141] Step 5-1: Pyramid Encoder Training: Image Stack The input is a pyramid encoder, which extracts multi-scale image features using a convolutional pyramid structure. The encoder employs... Layered convolutional structures are used to extract multi-scale feature representations of images:

[0142]

[0143] in, Representing the Multi-scale feature maps output after convolutional layers. The input image is used as the feature map of layer 0 of the network. No. The layer contains convolutional and batch-normalized convolutional operations. The LeakyReLU activation function (negative slope) =0.2), This represents the index of the convolutional layer used in the network. The final high-dimensional feature is represented as:

[0144]

[0145] in, It is the last layer of the network (the first layer) The output feature map generated by the layer, To output the number of feature channels, This represents the spatial dimensions (height and width) of the output feature map. This represents the high-dimensional, multi-scale image features ultimately extracted by the entire encoder.

[0146] Step 5-2: VAE encoder training: training high-dimensional features Parameters mapped to the latent distribution through fully connected layers:

[0147] (1) Gaussian distribution approximation: The VAE encoder learns the mean and log-variance parameters through a neural network:

[0148]

[0149]

[0150] in, and These are potential embedded variables. The mean and variance, Indicates that the encoder processes the input. The hidden representation obtained afterwards and For different fully connected layers. The posterior distribution defined by the encoder is:

[0151]

[0152] in, It is a latent embedded variable The conditional distribution, given the input The subsequent distribution.

[0153] (2) Reparameterization technique: To make the sampling process differentiable, reparameterization technique is used to generate latent embedded variables. :

[0154]

[0155] in, , This represents element-wise multiplication. Ultimately, the latent embedded variable... It is further reshaped into a feature map format for use in subsequent decoding processes.

[0156] Step 5-3: Fusion Layer Network Training: The goal of this module is to fuse multiple latent variables into a single latent representation. Nonlinear fusion of multiple latent variables is achieved through a densely connected network structure. The specific process is as follows:

[0157] (1) Multiple latent variables concatenation: ... The input latent embedded variables are concatenated along the channel dimension:

[0158]

[0159] in, It involves multiple potential embedded variables. The result of concatenating along the channel dimension forms a new latent variable. Each It is a size of The latent variables, the concatenated result Dimensions ,in , representing the total number of channels.

[0160] (2) Densely Connected Converged Network: Adopting Feature fusion is performed using a densely connected layer structure.

[0161]

[0162]

[0163] in, This indicates the total number of network layers, with a value of 5. Represents the index of the current layer; The input feature map is composed of multiple features concatenated together. Convolution operation in layer l; LeakReLU activation function (negative slope) =0.2); For the first The output feature map of the layer. Finally, the fused latent variable features are:

[0164]

[0165] in, Representing the The output feature map of the layer, It is the final fusion feature, which is composed of the outputs of all layers.

[0166] (3) Latent variable features are compressed into a single latent variable: The fused features are compressed into the target dimension through multi-layer convolution:

[0167]

[0168]

[0169] in, To compress the total number of layers in the network. Represents the index of the current compression layer. ; The input features are the fused features output from the previous stage; No. Layers are used for convolutional operations to compress features; The activation function is LeakyReLU (negative slope). =0.2); After the first Intermediate features obtained after layer compression. Finally, a single latent variable is obtained after fusion:

[0170]

[0171] in, The final output, a compressed single latent variable, is the output of the last layer. .

[0172] Step 5-4: U-Net Decoder: Reconstructs the original image from the fused latent variables. The U-Net architecture is used to implement the mapping from the latent space to the image space. The specific process is as follows:

[0173]

[0174] in, The U-Net decoder performs upsampling step by step through a symmetrical encoder-decoder path and skip connections, ultimately outputting a reconstructed image with pixel values ​​in the range [0, 1].

[0175] Step 5-5: Loss Function: The training objective is to minimize the weighted sum of the reconstruction loss and the KL divergence loss. To measure the difference between the decoded result and the original input image, the reconstruction loss is composed of both structural similarity loss and mean squared error loss. Specifically, the structural similarity loss... Defined as:

[0176]

[0177] in, Represents the original input image. This represents the image reconstructed by the decoder. This represents the structural similarity index, used to evaluate the consistency between the reconstructed image and the original image in terms of structure, brightness, and contrast. Mean squared error loss. Defined as:

[0178]

[0179] in, This represents the Frobenius norm, used to measure the reconstruction error of an image at the pixel level. , , , These represent the number of images to be fused, the number of channels in each image, and their height and width, respectively. Normalization is used to ensure the loss value remains numerically stable across different image pixel distributions. The total reconstruction loss is a weighted sum of the two values:

[0180]

[0181] Among them, weight The numerical scale is used to balance the two losses. Simultaneously, to avoid disordered latent spatial distribution and model overfitting, a KL divergence term is introduced into the objective function to measure the approximate posterior. With standard Gaussian prior Differences:

[0182]

[0183] in, and They represent the latent variables respectively. The mean and standard deviation of each dimension The KL divergence, representing the dimension of the latent space, can help the latent representation approximate a standard normal distribution, thereby improving its distribution regularity and separability. The final overall optimization objective function... The weighted sum of the two:

[0184]

[0185] in, This is a balancing coefficient used to adjust the relative importance of the reconstruction loss and the KL divergence loss. When... When the value is close to 0, the model prioritizes image reconstruction quality; when... When the value approaches 1, the model places more emphasis on the regularization constraints of the latent space.

[0186] During training, the Adam optimizer is used to jointly update the model parameters. In each iteration, forward propagation is performed to calculate the loss mentioned above, and backpropagation is used to update the parameters, thereby enabling the model to learn an effective encoder-fusion-decoder mapping and achieve high-quality multi-image fusion.

[0187] Step 6: Image Fusion Output: After model training is complete, the reconstructed images generated by the U-Net decoder will be fused and output. As the final output.

[0188] Figure 6 This paper demonstrates the effectiveness of a multi-frame photoacoustic image reconstruction method based on optical flow alignment and depth feature fusion in simulation experiments. Simulation results show that the method can effectively recover image details and reduce distortion caused by motion and noise under ideal conditions. This further verifies that the framework has good detail recovery capability and robustness in dynamic photoacoustic imaging, providing a theoretical basis and technical support for subsequent practical applications.

Claims

1. A multi-frame photoacoustic image reconstruction method based on optical flow alignment and depth feature fusion, characterized by The method comprises the following steps: Step S1, photoacoustic signal data acquisition; Step S2, photoacoustic tomographic image reconstruction; Step S3, optical flow calculation and image alignment: Step S3-1, calculation of optical flow field: select a frame from the initial image sequence as the reference image, and calculate the optical flow field between the remaining target images and the reference image based on the optical flow method, the specific steps are as follows: Input reference image and target image Each pixel neighborhood is modeled using a quadratic polynomial: A polynomial expansion is performed, fitting a quadratic polynomial to each pixel point in both images, in the form: ； wherein is a coordinate in local coordinates, is a symmetric matrix, is a vector, is a scalar, denotes the transpose operation; On the constructed image pyramid, a coarse-to-fine estimation is performed, starting from the coarsest level, an initial displacement field is assumed, at each level, based on the principle that the polynomial coefficients should remain unchanged under motion, the displacement vector of each pixel is solved by solving a weighted least squares system constructed by all pixels in the neighborhood The process satisfies the brightness constancy constraint: ； wherein, is the pixel intensity of the reference image in global space coordinates and time instant is the pixel intensity of the target image in the corresponding shifted position, is the pixel intensity of the target image in the corresponding shifted position, is the global space coordinate of the pixel in the image, is the time of the current image frame, are the horizontal and vertical displacement of the pixel between the two frames, respectively. Through iterative optimization, the displacement field is continuously updated to minimize the intensity difference between the two images, and the objective function is: ； wherein, the reference image at the pixel has a luminance value; the target image at the shifted position has a luminance value; the index of the representative pixel; the displacement that minimizes the error ; The hierarchical estimated optical flow is up-sampled to a preset number of levels as an initial value, and the above optimization process is repeated until the original resolution is reached, and finally, a dense optical flow field is output ; Step S3-2, image alignment: the calculated optical flow field is applied to each target image to perform geometric transformation, thereby obtaining an image sequence that is strictly aligned in space with the reference image; Step S4, depth feature fusion network training: Step S4-1, data preparation: the aligned target image and the reference image are merged in the channel and normalized to form a network input; Step S4-2, pyramid encoder training; Step S4-3, VAE encoder training; Step S4-4, fusion layer network training: the fusion layer network receives all latent embedding variables from the VAE encoder The information is integrated by the convolution layer with dense connection and reflection padding structure, and a unified and information-rich fusion latent representation is output. The specific steps are as follows: Multi-potential variable concatenation: concatenating one input potential embedding variable along the channel dimension: ； wherein, is the result of concatenating the plurality of potential embedding variables along the channel dimension, denotes the sum of the number of channels, represents the spatial size of the output feature map; Dense connection fusion network: adopt Layer dense connection structure for feature fusion: ； ； wherein, represents the total number of network layers, and takes a value of 5, represents the index of the current layer, is the The layer contains a convolution operation of convolution and batch normalization, is a LeakyReLU activation function, is the The output feature map of the layer, finally, the latent variable feature is fused as: ； wherein is the final fused feature; Latent variable feature compression into a single latent variable: the fusion feature is compressed into the target dimension through multiple layers of convolution: ； ； wherein, is the total number of layers of the compressed network, represents the index of the current compressed layer, the first convolutional operation for feature compression, the first intermediate features obtained after layer compression, and finally a fused single latent variable: ； wherein, the final output, compressed single latent variable, i.e. the output of the last layer ; Step S4-5, U-Net decoder training; Step S4-6, loss function calculation: loss function consists of three parts: structural similarity loss, mean square error loss and KL divergence loss Step S4-7, optimization process: using the Adam optimizer, the parameters of the pyramid encoder, the VAE encoder, the fusion layer network and the U-Net decoder are jointly optimized by minimizing the loss function.

2. The multi-frame photoacoustic image reconstruction method based on light flow alignment and depth feature fusion according to claim 1, characterized in that The specific steps of step 2 are as follows: Step S2-1, setting imaging parameters: setting the radius of the annular ultrasonic probe, the number of ultrasonic transducers, the detection angle, the size of the imaging physical region, the size of the pixel grid, the sampling frequency of the data acquisition device, the sound speed of sound wave propagation and the aperture size parameters in the imaging process; Step S2-2, image reconstruction using DAS imaging method: applying the DAS imaging method to reconstruct the image sequence.

3. The multi-frame photoacoustic image reconstruction method based on fusion of optical flow alignment and depth features according to claim 2, characterized in that In the DAS imaging method in step S2-2, the relationship between the beamforming signal and the received signal is defined as: wherein is the number of ultrasound transducers, denotes a beamformed signal at a position , denotes a signal received by the th transducer element at time , denotes the time it takes for a sound wave originating from the position to travel to the th transducer element, and is given by ； wherein representing a point distance between the first transducer elements, is the speed of sound.

4. The multi-frame photoacoustic image reconstruction method based on fusion of optical flow alignment and depth features according to claim 1, characterized in that The specific steps of step 3-2 are as follows: Obtaining an optical flow field After that, by geometrically deforming the target image to make it spatially align with the reference image , for each coordinate position in the aligned image , its corresponding sub-pixel coordinate in the source image is found according to the optical flow field , the pixel intensity value at this position is calculated from using bilinear interpolation and is assigned to , after traversing all pixels, the aligned image that is exactly registered with the reference image is finally generated.

5. The multi-frame photoacoustic image reconstruction method based on fusion of optical flow alignment and depth features according to claim 1, characterized in that In step 4-2, the processed image is input into the pyramid encoder to extract multi-scale features, and the high-quality feature map is obtained, and the encoder parameters are optimized through gradient back propagation, and the specific steps are as follows: stack of images input pyramid encoder, which extracts multi-scale image features using a convolution pyramid structure, the encoder adopts layer convolution structure to extract multi-scale feature representation of the image: ； wherein, represent the first multi-scale feature maps output after the convolution of the layer, is the input image as the 0th layer feature map of the network, the first contains a convolution operation with batch normalization, is a LeakyReLU activation function, represents the index of the convolution layer used in the network, and the final high-dimensional feature representation is: ； wherein, is the network first layer generated output feature map, is the output feature channel number, represents the spatial size of the output feature map, represents the high-dimensional multi-scale image features finally extracted by the entire encoder.

6. The multi-frame photoacoustic image reconstruction method based on fusion of optical flow alignment and depth features according to claim 5, characterized in that In step 4-3, the feature map obtained in step S4-2 is input into the encoder of the variational autoencoder VAE, the posterior distribution of the latent space is modeled using Gaussian distribution, the mean vector and the standard deviation vector of each feature are learned, and the reparameterization technique is used to realize the differentiability of the sampling process, which is convenient for gradient optimization, and the specific steps are as follows: Gaussian distribution approximation: the VAE encoder learns the mean and log-variance parameters through a neural network: ； ； where, and are the mean and variance of the latent embedding variable denotes the hidden representation obtained after the encoder processes the input and are different fully connected layers; the posterior distribution defined by the encoder is:​​ ； where, is the conditional distribution of the latent embedding variable given the input the distribution after the embedding Reparameterization trick: To make the sampling process differentiable, the reparameterization technique is used to generate latent variables : ； wherein , denotes an element-wise multiplication.

7. The multi-frame photoacoustic image reconstruction method based on fusion of optical flow alignment and depth features of claim 1, wherein In step 4-5, the decoder receives the fused latent variable, gradually reconstructs the spatial structure of the feature through convolution, upsampling and skip connection, enhances the non-linear expression ability by using the ReLU activation function, and finally outputs the fused high-quality image, and the specific steps are as follows: The fused latent variable is reconstructed into the original image, and the U-Net architecture is used to realize the mapping from the latent space to the image space, and the specific process is as follows: ； wherein, is the U-Net decoder, which gradually upsamples through a symmetric encoding-decoding path and skip connections, finally outputting a reconstructed image with pixel values in the range [0, 1].

8. The multi-frame photoacoustic image reconstruction method based on fusion of optical flow alignment and depth features according to claim 7, characterized in that In steps 4-6, the structural similarity loss is defined as: ； wherein, denotes the original input image, denotes the decoder reconstructed image, denotes the structural similarity index; mean squared error loss is defined as: ； wherein, denotes the Frobenius norm, , , , denote the number of images to be fused, the number of channels of the images, the height and the width, respectively. KL divergence loss is defined as: ； where and denote the mean and standard deviation of the latent variable dimension, respectively, denotes the dimension of the latent space. Loss function is defined as: ； ； wherein, is a balancing coefficient, is the total reconstruction loss, is a weight.

Images