A photoacoustic remote sensing microscopy imaging method and imaging system based on parallel processing acceleration

Abstract

The present invention is applicable to the technical field of optical imaging, and provides a photoacoustic remote sensing microscopic imaging method and an imaging system based on parallel processing acceleration. A graphics processing unit (GPU) is used to introduce a unified device architecture based on parallel processing into photoacoustic remote sensing microscopic imaging to improve data processing speed. The data transmission form of data stream is used to break through the limitation of the traditional PARS signal acquisition method subject to the reloading time of the acquisition card. A multi-atomic lock algorithm is used to avoid memory competition caused by multiple threads simultaneously accessing a shared memory location in concurrent threads, thereby greatly improving efficiency.

Description

Technical Field

[0001] The invention belongs to the technical field of optical imaging, in particular to an acceleration method and system for photoacoustic remote sensing microscopic imaging. Background Art

[0002] Photoacoustic remote sensing (PARS) technology uses low-coherent light to perform non-contact photoacoustic detection under the principle of incoherence, which can avoid phase noise in principle, and has extremely high resolution and higher signal-to-noise ratio than traditional contact photoacoustic microscopy. The miniaturization, non-contact, label-free and fast imaging characteristics of PARS make this technology particularly suitable for application scenarios such as intraoperative fast microscopic imaging and fundus vascular functional imaging.

[0003] The fast imaging characteristics of PARS place higher speed requirements on the system's laser repetition rate, signal acquisition, and data processing. The higher the laser repetition rate, the more PARS signals are excited per unit time, the larger the amount of data collected, and the longer it takes to process the data. When the laser repetition rate is greater than 200KHz, the theoretical imaging frame rate will be greater than two frames per second. A microscopic image with a 100-micron imaging field of view will need to generate several GB of data, and the data acquisition and processing must be completed in a time far less than a second.

[0004] The traditional PARS signal acquisition method is triggered acquisition, that is, the data acquisition card triggers an acquisition and transmission action under the specified trigger conditions, and the acquired data is non-continuous and strictly related to the trigger sequence. After each trigger, the acquisition card needs to be reloaded to welcome the next trigger. In this acquisition method, when the interval time of the trigger sequence is less than the reloading time of the acquisition card, the acquired data will be missed and delayed. When processing data, traditional PARS systems usually use CPU calculations, but the CPU serial processing of large amounts of data (several GB level) is time-consuming, which limits the imaging frame rate of the PARS system. Summary of the invention

[0005] In order to overcome the shortcomings and deficiencies of the prior art, the present invention designs a photoacoustic remote sensing microscopy imaging method and an imaging system based on parallel processing acceleration based on GPU parallel processing technology.

[0006] The photoacoustic remote sensing microscopic imaging method and imaging system based on parallel processing acceleration include: a data acquisition module, a data extraction module, a data preprocessing module, an extreme value finding module, an interpolation module and an imaging module.

[0007] In the data acquisition module, in the data stream mode, the acquisition card acquires PARS signals, trigger signals, and position information of the fast axis and the slow axis, and the data of the four acquisition channels are in one data stream; the CPU opens up a first buffer A and a second buffer B in the memory, and the space size of the buffers is 4*fs*t, wherein fs represents the sampling frequency of the acquisition card, and t is the sampling time; the first buffer A is set as a receiving area, and the second buffer B is set as a processing area; the receiving area receives the data flowing out from the acquisition card, and the processing area (if not empty) copies the data in the area to the GPU video memory, and the allocated video memory space size is 4*fs*t; if the receiving area is full and the processing area is processed, the receiving area and the processing area are switched to continue receiving and processing data.

[0008] During data acquisition, in data stream mode, there is only a single trigger in data acquisition through the sampling frequency of the acquisition card. The acquisition card hardware does not need to be reloaded and is no longer limited by the reload time.

[0009] In the data extraction module, a two-dimensional thread is allocated in the GPU, the thread size is 1024, and the block size is 4*fs*t / 1024, which is used to perform effective data extraction on the sampled data.

[0010] Since the data acquisition module continuously acquires signals according to the sampling frequency of the acquisition card, the sampled data acquired by the acquisition card includes the excited and unexcited PARS signals. Therefore, both the excited and unexcited PARS signals can be completely acquired without the problem of data loss. The continuous acquisition of PARS signals and the triggered acquisition of PARS signals are as follows: Figure 4 However, since the unexcited PARS signal is collected, it is necessary to perform effective data extraction on the sampled data. The specific extraction method is as follows:

[0011] In the data stream, the index values ​​of the acquisition channel data at the same time are adjacent. The data of the four acquisition channels are divided with four data as a cycle; the data in the trigger signal data whose data value is less than the threshold is set to zero. Here, the threshold is set to 0.3. Other thresholds can also be set according to actual needs, as follows:

[0012]

[0013] Calculate the rising edge position of the processed data, and the judgment formula is:

[0014]

[0015] Where f(n) is the original trigger signal data value, f′(n) is the processed trigger signal data value, n is the index value, and the index value that meets the above conditions is the position where the rising edge occurs.

[0016] In the data preprocessing module, the valid data is converted into corresponding position data and PARS signal peak-to-peak data, and only the valid data in the sampled data is calculated. The specific calculation method is as follows:

[0017] Since the two-dimensional galvanometer realizes plane scanning by driving the reflective lens deflection through the XY axis motor, it is generally used for smaller scanning areas, and the scanning frequency is much smaller than the sampling frequency. After a single trigger, the subsequent m valid data can be considered as signals at the same position. The average value of the m valid data is calculated as the signal value at the time of a single trigger to calculate the position data. The calculation formula is:

[0018]

[0019] Where X and Y represent the position information of the fast axis and the slow axis respectively, d is the focal length of the objective lens, k is the voltage to radian conversion ratio, sum x and sum y are the sum of m valid data of the fast axis and the slow axis, respectively, x and z y are the initial position voltages of the fast axis and the slow axis, respectively.

[0020] The peak-to-peak value is obtained by calculating the difference between the maximum and minimum values ​​of the subsequent m valid data of the PARS signal. The formula is:

[0021] V PP =V max -V min (4)

[0022] Where V PP is the peak-to-peak value, V max is the maximum value of the subsequent m valid data, V min It is the minimum value of the subsequent m valid data.

[0023] In the extreme value seeking module, a reduction algorithm based on an alternating strategy is used to obtain the extreme value of the position data, a two-dimensional thread is allocated in the GPU, and the product of the thread size and the block size is used as the alternating factor. The GPU is used to calculate the position data to obtain an intermediate result array, and then the result array is copied to the CPU; in the CPU, the final extreme value of the position data is calculated by serial processing.

[0024] In the gridding module, the position information of the X-axis and the Y-axis needs to be scaled to the total pixel width of the final image required, and the calculation formula is:

[0025]

[0026] where x i and iare the X and Y coordinate values ​​of the i-th scanning point, respectively. max 、x min ,y max and min are the maximum and minimum values ​​of the X and Y coordinates obtained in the extreme value step. x and N y are the number of grids that the image is divided into in the X and Y directions (the final image pixel width). xi and n yi The grid coordinates corresponding to the i-th scanning point in the X and Y directions respectively; according to the calculated coordinate information, a multiple atomic lock algorithm is used to apply for each grid's own atomic lock to synchronize the threads in the thread block. When multiple threads apply to access the same address, the atomic lock limits access to only one thread at a time, and the remaining threads wait for the atomic lock to be released; in the two-dimensional thread, GPU statistics are used to obtain the number of scanning points contained in each grid, and the peak-to-peak value data in the valid data is used as the signal strength value of each scanning point. After each thread accesses, the result is visible to all threads.

[0027] In the interpolation module, for a grid containing multiple scanning points, the average signal strength value is calculated according to the number of scanning points as the signal strength value of the grid; for a grid not containing scanning points, a linear interpolation method is used to fill it, and the calculation formula is:

[0028]

[0029] Where i and j are the row and column coordinates of the grid, respectively. i+1,j 、V i,j+1 、V i-1,j 、V i,j-1 is the sum of the signal strength values ​​of each grid, N i+1,j 、N i,j+1 、N i-1,j 、N i,j-1 is the number of scan points of the respective grid, v i,j is the average signal value of the grid. The signal strength value data of all grid points calculated are used as grid data.

[0030] In the imaging module, the grid data obtained after processing is used to find the extreme value using a reduction algorithm, and after grayscale normalization processing, the calculation formula is:

[0031]

[0032] Where i and j are the grid row and column coordinates, I i,j is the pixel value after grayscale normalization, v i,j is the grid signal value, v max 、v min are the maximum and minimum signal strength values ​​of all grids.

[0033] The present invention also includes a photoacoustic remote sensing microscopic imaging method and an imaging method based on parallel processing acceleration.

[0034] Compared with the prior art, the present invention has the following advantages and beneficial effects:

[0035] 1. The present invention is based on a photoacoustic remote sensing microscopy imaging method and system accelerated by parallel processing. The unified device architecture based on parallel processing is introduced into the photoacoustic remote sensing microscopy imaging by using a graphics processing unit (GPU). The GPU is more adept at parallel computing of large amounts of data and floating-point data. The problem that as the imaging speed increases in photoacoustic remote sensing microscopy, the data volume increases sharply, which in turn leads to an increase in data processing difficulty and processing time is solved. Compared with the data imaging speed of serial processing, the imaging speed of parallel processing is increased by dozens of times. Data collection and processing are performed simultaneously, and the data processing time is much shorter than the collection time, so that the imaging speed is only limited by the data collection time (determined by the repetition rate of the pulsed laser), and the actual imaging speed is closer to the theoretical speed.

[0036] 2. The present invention adopts the data transmission form of data stream, breaking through the limitation of the traditional PARS signal acquisition method subject to the reloading time of the acquisition card. Since the PARS signal length is very short and the trigger frequency is very high, as the imaging speed increases, the effective data in the data stream increases, and the data processing efficiency is higher. The acceleration effect of the present invention is more obvious, which is a more applicable PARS signal transmission and processing method, and meets the development requirements of photoacoustic remote sensing microscopy imaging technology in rapid imaging.

[0037] 3. Since the present invention adopts a photoacoustic remote sensing microscopy imaging method based on parallel processing acceleration, all complex data processing and image processing processes are completed in the GPU. Since the parallel processing acceleration method can be run under relatively low hardware configuration, it does not require the use of high-end hardware equipment, thereby reducing the cost of the system.

[0038] 4. The present invention adopts a multiple atomic lock algorithm to avoid memory contention caused by multiple threads accessing a shared memory location at the same time in concurrent threads, which greatly improves efficiency and ensures the atomicity of the operation, that is, the execution process of "read-modify-write" a variable is the smallest operation unit and cannot be further decomposed. Therefore, during the execution of atomic operations, other parallel threads cannot read and write the variable, thereby ensuring data consistency and correctness, and also improving concurrency performance. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] In order to more clearly illustrate the technical solutions in the embodiments of the present invention, the drawings required for use in the embodiments or the description of the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying creative labor.

[0040] Figure 1 A flow chart of the photoacoustic remote sensing microscopy imaging method based on parallel processing acceleration of the present invention;

[0041] Figure 2 This is a block diagram of the photoacoustic remote sensing microscopy imaging system based on parallel processing acceleration;

[0042] Figure 3 It is the calculation flow chart of the data processing module;

[0043] Figure 4 It is the timing diagram of data collection;

[0044] Figure 5 The figure is a comparison diagram of the imaging results of the imaging method of the present invention and the imaging method without acceleration. DETAILED DESCRIPTION

[0045] In order to enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be clearly and completely described below in conjunction with the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all of the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without making creative work are within the scope of protection of the present application.

[0046] Reference to "embodiments" in this application means that a particular feature, structure, or characteristic described in conjunction with the embodiments may be included in at least one embodiment of the present application. The appearance of the phrase in various locations in the specification does not necessarily refer to the same embodiment, nor is it an independent or alternative embodiment that is mutually exclusive with other embodiments. It is explicitly and implicitly understood by those skilled in the art that the embodiments described in this application may be combined with other embodiments.

[0047] like Figure 1-2 As shown, the photoacoustic remote sensing microscopy imaging method and imaging system based on parallel processing acceleration include: a data acquisition module, a data extraction module, a data preprocessing module, an extreme value finding module, an interpolation module and an imaging module.

[0048] In the data acquisition module, in the data stream mode, the acquisition card acquires PARS signals, trigger signals, and position information of the fast axis and the slow axis, and the data of the four acquisition channels are in one data stream; the CPU opens up a first buffer A and a second buffer B in the memory, and the space size of the buffers is 4*fs*t, wherein fs represents the sampling frequency of the acquisition card, and t is the sampling time; the first buffer A is set as a receiving area, and the first buffer A is set as a processing area; the receiving area receives the data flowing out from the acquisition card, and the processing area (if not empty) copies the data in the area to the GPU video memory, and the allocated video memory space size is 4*fs*t; if the receiving area is full and the processing area is processed, the receiving area and the processing area are switched to continue receiving and processing data.

[0049] During data acquisition, in data stream mode, there is only a single trigger in data acquisition through the sampling frequency of the acquisition card. The acquisition card hardware does not need to be reloaded and is no longer limited by the reload time.

[0050] In the data extraction module, a two-dimensional thread is allocated in the GPU, the thread size is 1024, and the block size is 4*fs*t / 1024, which is used to perform effective data extraction on the sampled data.

[0051] Since the data acquisition module continuously acquires signals according to the sampling frequency of the acquisition card, the sampled data acquired by the acquisition card includes the excited and unexcited PARS signals. Therefore, both the excited and unexcited PARS signals can be completely acquired without the problem of data loss. The continuous acquisition of PARS signals and the triggered acquisition of PARS signals are as follows: Figure 4 However, since the unexcited PARS signal is collected, it is necessary to perform effective data extraction on the sampled data. The specific extraction method is as follows:

[0052] In the data stream, the index values ​​of the acquisition channel data at the same time are adjacent. The data of the four acquisition channels are divided with four data as a cycle; the data in the trigger signal data that is less than the threshold is set to zero. Here, the threshold is set to 0.3. Other thresholds can also be set according to actual needs, as follows:

[0053]

[0054] Calculate the rising edge position of the processed data, and the judgment formula is:

[0055]

[0056] Where f′(n) is the data value of the trigger signal, n is the index value, and the index value that meets the above conditions is the position where the rising edge occurs.

[0057] In the data preprocessing module, the valid data is converted into corresponding position data and PARS signal peak-to-peak data, and only the valid data in the sampled data is calculated. The specific calculation method is as follows:

[0058] Since the two-dimensional galvanometer realizes plane scanning by driving the reflective lens deflection through the XY axis motor, it is generally used for smaller scanning areas, and the scanning frequency is much smaller than the sampling frequency. After a single trigger, the subsequent m valid data can be considered as signals at the same position. The average value of the m valid data is calculated as the signal value at the time of a single trigger to calculate the position data. The calculation formula is:

[0059]

[0060] Where X and Y represent the position information of the fast axis and the slow axis respectively, d is the focal length of the objective lens, k is the voltage to radian conversion ratio, sum x and sum y are the sum of m valid data of the fast axis and the slow axis, respectively, x and z y are the initial position voltages of the fast axis and the slow axis, respectively.

[0061] The peak-to-peak value is obtained by calculating the difference between the maximum and minimum values ​​of the subsequent m valid data of the PARS signal. The formula is:

[0062] V PP =V max -V min (4)

[0063] Where V PP is the peak-to-peak value, V max is the maximum value of the subsequent m valid data, V min It is the minimum value of the subsequent m valid data.

[0064] In the extreme value seeking module, a reduction algorithm based on an alternating strategy is used to obtain the extreme value of the position data, a two-dimensional thread is allocated in the GPU, and the product of the thread size and the block size is used as the alternating factor. The GPU is used to calculate the position data to obtain an intermediate result array, and then the result array is copied to the CPU; in the CPU, the final extreme value of the position data is calculated by serial processing.

[0065] In the gridding module, the position information of the X-axis and the Y-axis needs to be scaled to the total pixel width of the final image required, and the calculation formula is:

[0066]

[0067] where x i and iare the X and Y coordinate values ​​of the i-th scanning point, respectively. max 、x min ,y max and min are the maximum and minimum values ​​of the X and Y coordinates obtained in the extreme value step. x and N y are the number of grids that the image is divided into in the X and Y directions (the final image pixel width). xi and n yi The grid coordinates corresponding to the i-th scanning point in the X and Y directions respectively; according to the calculated coordinate information, use the multiple atomic lock algorithm to apply for each grid's own atomic lock, synchronize the threads in the thread block, and when multiple threads apply to access the same address, the atomic lock limits access to only one thread at a time, and the remaining threads wait for the atomic lock to be released; in the two-dimensional thread, use GPU statistics to obtain the number of scanning points contained in each grid, and use the peak-to-peak value data in the valid data as the signal strength value of each scanning point. After each thread accesses, the result will be visible to all threads.

[0068] In the interpolation module, for a grid containing multiple scanning points, the average signal strength value is calculated according to the number of scanning points as the signal strength value of the grid; for a grid not containing scanning points, a linear interpolation method is used to fill it, and the calculation formula is:

[0069]

[0070] Where i and j are the row and column coordinates of the grid, respectively. i+1,j 、V i,j+1 、V i-1,j 、V i,j-1 is the sum of the signal strength values ​​of each grid, N i+1,j 、N i,j+1 、N i-1,j 、N i,j-1 is the number of scan points of the respective grid, v i,j is the average signal value of the grid. The signal strength value data of all grid points calculated are used as grid data.

[0071] In the imaging module, the grid data obtained after processing is used to find the extreme value using a reduction algorithm, and after grayscale normalization processing, the calculation formula is:

[0072]

[0073] Where i and j are the grid row and column coordinates, I i,j is the pixel value after grayscale normalization, v i,j is the grid signal value, v max 、v minare the maximum and minimum signal strength values ​​of all grids; the results are then displayed.

[0074] like Figure 3 As shown, the photoacoustic remote sensing microscopy imaging method based on parallel processing acceleration includes the following steps:

[0075] Step 1: In data stream mode, the acquisition card acquires PARS signals, trigger signals, and the position information of the fast axis and slow axis. The data of the four acquisition channels are in one data stream. The CPU opens up a first buffer A and a second buffer B in the memory. The space size of the buffers is 4*fs*t, where fs represents the sampling frequency of the acquisition card and t is the sampling time. The first buffer A is set as a receiving area and the second buffer A is set as a processing area. The receiving area receives the data flowing out of the acquisition card, and the processing area (if not empty) copies the data in this area to the GPU video memory, and allocates a video memory space size of 4*fs*t. If the receiving area is full and the processing area has finished processing, the receiving area and the processing area are switched to continue receiving and processing data.

[0076] During data acquisition, in data stream mode, there is only a single trigger in data acquisition through the sampling frequency of the acquisition card. The acquisition card hardware does not need to be reloaded and is no longer limited by the reload time.

[0077] Step 2: Allocate a two-dimensional thread in the GPU with a thread size of 1024 and a block size of 4*fs*t / 1024 for effective data extraction of the sampled data.

[0078] Since the data acquisition module continuously acquires signals according to the sampling frequency of the acquisition card, the sampled data acquired by the acquisition card includes the excited and unexcited PARS signals. Therefore, both the excited and unexcited PARS signals can be completely acquired without the problem of data loss. The continuous acquisition of PARS signals and the triggered acquisition of PARS signals are as follows: Figure 4 However, since the unexcited PARS signal is collected, it is necessary to perform effective data extraction on the sampled data. The specific extraction method is as follows:

[0079] In the data stream, the index values ​​of the acquisition channel data at the same time are adjacent. The data of the four acquisition channels are divided with four data as a cycle; the data in the trigger signal data whose data value is less than the threshold is set to zero. Here, the threshold is set to 0.3. Other thresholds can also be set according to actual needs, as follows:

[0080]

[0081] Calculate the rising edge position of the processed data, and the judgment formula is:

[0082]

[0083] Where f(n) is the original trigger signal data value, f′(n) is the processed trigger signal data value, n is the index value, and the index value that meets the above conditions is the position where the rising edge occurs.

[0084] Step 3: Convert the valid data into the corresponding position data and PARS signal peak-to-peak data, and only calculate the valid data in the sampled data. The specific calculation method is as follows:

[0085] Since the two-dimensional galvanometer realizes plane scanning by driving the reflective lens deflection through the XY axis motor, it is generally used for smaller scanning areas, and the scanning frequency is much smaller than the sampling frequency. After a single trigger, the subsequent m valid data can be considered as signals at the same position. The average value of the m valid data is calculated as the signal value at the time of a single trigger to calculate the position data. The calculation formula is:

[0086]

[0087] Where X and Y represent the position information of the fast axis and the slow axis respectively, d is the focal length of the objective lens, k is the voltage to radian conversion ratio, sum x and sum y are the sum of m valid data of the fast axis and the slow axis, respectively, x and z y are the initial position voltages of the fast and slow axes, respectively;

[0088] The peak-to-peak value is obtained by calculating the difference between the maximum and minimum values ​​of the subsequent m valid data of the PARS signal. The formula is:

[0089] V PP =V max -V min (4)

[0090] Where V PP is the peak-to-peak value, V max is the maximum value of the subsequent m valid data, V min It is the minimum value of the subsequent m valid data;

[0091] Step 4: In the extreme value seeking module, a reduction algorithm based on an alternating strategy is used to obtain the extreme value of the position data, a two-dimensional thread is allocated in the GPU, and the product of the thread size and the block size is used as the alternating factor. The GPU is used to calculate the position data to obtain an intermediate result array, and then the result array is copied to the CPU; in the CPU, the final extreme value of the position data is calculated by serial processing.

[0092] Step 5: Grid the X-axis and Y-axis position information to the required total image pixel width. The calculation formula is:

[0093]

[0094] where x i and i are the X and Y coordinate values ​​of the i-th scanning point, respectively. max 、x min ,y max and min are the maximum and minimum values ​​of the X and Y coordinates obtained in the extreme value step. x and N y are the number of grids that the image is divided into in the X and Y directions (the final image pixel width). xi and n yi The grid coordinates corresponding to the i-th scanning point in the X and Y directions respectively; according to the calculated coordinate information, a multiple atomic lock algorithm is used to apply for each grid's own atomic lock to synchronize the threads in the thread block. When multiple threads apply to access the same address, the atomic lock limits access to only one thread at a time, and the remaining threads wait for the atomic lock to be released; in the two-dimensional thread, GPU statistics are used to obtain the number of scanning points contained in each grid, and the peak-to-peak value data in the valid data is used as the signal strength value of each scanning point. After each thread accesses, the result is visible to all threads.

[0095] Step 6: Interpolate the grid. For grids containing multiple scanning points, calculate the average signal strength value according to the number of scanning points as the signal strength value of the grid; for grids that do not contain scanning points, use the linear interpolation method to fill them. The calculation formula is:

[0096]

[0097] Where i and j are the row and column coordinates of the grid, respectively. i+1,j 、V i,j+1 、V i-1,j 、V i,j-1 is the sum of the signal strength values ​​of each grid, N i+1,j 、N i,j+1 、N i-1,j 、N i,j-1 is the number of scan points of the respective grid, v i,j is the average signal value of the grid. The signal strength value data of all grid points calculated are used as grid data.

[0098] Step 7: Use the grid data obtained in step 6 to find the extreme value using the reduction algorithm, and obtain the final image after grayscale normalization. The calculation formula is:

[0099]

[0100] Where i and j are the grid row and column coordinates, I i,j is the pixel value after grayscale normalization, v i,j is the grid signal value, v max 、v min are the maximum and minimum signal strength values ​​of all grids.

[0101] The imaging effects after and without the accelerated parallel processing proposed in this scheme are shown in the figure below. Figure 5 As shown. Figure 5 It can be seen that under the continuous acquisition method, the reconstruction algorithms based on CPU serial processing and GPU parallel processing obtained the same imaging results, and no signal loss occurred, but the accelerated reconstruction time was shortened by about 30 times.

[0102] It should also be noted that in this specification, terms such as "comprises", "includes" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device including a series of elements includes not only those elements, but also includes other elements not explicitly listed, or also includes elements inherent to such process, method, article or device. In the absence of more restrictions, an element defined by the sentence "comprises a ..." does not exclude the presence of other identical elements in the process, method, article or device including the element.

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

[0104] The above embodiments are preferred implementation modes of the present invention, but the implementation modes of the present invention are not limited to the above embodiments. Any other changes, modifications, substitutions, combinations, and simplifications that do not deviate from the spirit and principles of the present invention should be equivalent replacement methods and are included in the protection scope of the present invention.

[0105] The embodiments described above are only used to illustrate the technical solutions of the present invention, rather than to limit the same. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that the technical solutions described in the aforementioned embodiments may still be modified, or some of the technical features may be replaced by equivalents. Such modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention, and should all be included in the protection scope of the present invention.

Claims

1. A photoacoustic remote sensing microscopy imaging system based on parallel processing acceleration, characterized in that: The imaging system comprises: a data acquisition module, a data extraction module, a data preprocessing module, an extreme value finding module, a gridding module, an interpolation module and an imaging module; In the data acquisition module, in the data stream mode, the acquisition card acquires PARS signals, trigger signals, and position information of the fast axis and the slow axis, and the data of the four acquisition channels are in one data stream; the CPU opens up a first buffer A and a second buffer B in the memory; the first buffer A is set as a receiving area, and the first buffer B is set as a processing area; the receiving area receives the data flowing out from the acquisition card, and if the processing area is not empty, the data in the area is copied to the GPU video memory; if the receiving area is full and the processing area has finished processing, the receiving area and the processing area are switched to continue receiving and processing data; The data extraction module allocates two-dimensional threads in the GPU for performing effective data extraction on the sampled data; The index values ​​of the acquisition channel data at the same time are adjacent. The data of the four acquisition channels are divided with four data as a cycle. The data values ​​in the trigger signal data that are less than the threshold are set to zero. The threshold is set to 0.

3. The calculation formula is shown in (1): Calculate the rising edge position of the processed data, and the judgment formula is: Where f(n) is the original trigger signal data value, f′(n) is the processed trigger signal data value, n is the index value, and the index value that meets the above conditions is the position where the rising edge occurs; The data preprocessing module converts the valid data into corresponding position data and PARS signal peak-to-peak data, and only calculates the valid data in the sampled data. The specific calculation method is as follows: After a single trigger, the average value of m valid data is calculated as the signal value at the time of single trigger to calculate the position data. The calculation formula is: Where X and Y represent the position information of the fast axis and the slow axis respectively, d is the focal length of the objective lens, k is the voltage to radian conversion ratio, sum x and sum y are the sum of m valid data of the fast axis and the slow axis, respectively, x and z y are the initial position voltages of the fast and slow axes, respectively; The peak-to-peak value is obtained by calculating the difference between the maximum and minimum values ​​of the subsequent m valid data of the PARS signal. The formula is: V PP =V max -V min (4) Where V PP is the peak-to-peak value, V max is the maximum value of the subsequent m valid data, V min It is the minimum value of the subsequent m valid data; The extreme value seeking module obtains the extreme value of the position data based on the reduction algorithm of the alternating strategy, allocates two-dimensional threads in the GPU, uses the product of the thread size and the block size as the alternating factor, uses the GPU to calculate the position data to obtain an intermediate result array, and then copies the result array to the CPU; in the CPU, the final extreme value of the position data is calculated in a serial processing manner; The gridding module needs to scale the position information of the X-axis and Y-axis to the final required total image pixel width, and the calculation formula is: where x i and i are the X and Y coordinate values ​​of the i-th scanning point, respectively. max 、x min ,y max and min are the maximum and minimum values ​​of the X and Y coordinates obtained in the extreme value step; N x and N y are the number of grids that the image is divided into in the X and Y directions, i.e., the pixel width of the final image; n xi and n yi The grid coordinates corresponding to the i-th scanning point in the X and Y directions are respectively; according to the calculated coordinate information, a multi-atomic lock algorithm is used to apply for each grid's own atomic lock to synchronize the threads in the thread block. When multiple threads apply to access the same address, the atomic lock limits access to only one thread at a time, and the remaining threads wait for the atomic lock to be released; GPU statistics are used in the two-dimensional thread to obtain the number of scanning points contained in each grid, and the peak-to-peak value data in the valid data is used as the signal strength value of each scanning point. After each thread accesses, the result is visible to all threads; The interpolation module calculates the average signal strength value according to the number of scanning points as the signal strength value of the grid containing multiple scanning points; for the grid not containing scanning points, the linear interpolation method is used to fill it, and the calculation formula is: Where i and j are the row and column coordinates of the grid, respectively. i+1,j 、V i,j+1 、V i-1,j 、V i,j-1 is the sum of the signal strength values ​​of each grid, N i+1,j 、N i,j+1 、N i-1,j 、N i,j-1 is the number of scan points of the respective grid, v i,j is the average signal value of the grid; the signal strength value data of all grid points calculated are used as grid data; The imaging module uses a reduction algorithm to find the extreme value of the processed grid data, and performs grayscale normalization processing to obtain the final display image. The calculation formula is: Where i and j are the grid row and column coordinates, I i,j is the pixel value after grayscale normalization, v i,j is the grid signal value, v max 、v min are the maximum and minimum signal strength values ​​of all grids.

2. The imaging system as claimed in claim 1, wherein the space sizes of the first buffer A and the second buffer B are both 4*fs*t, wherein fs represents the sampling frequency of the acquisition card, and t represents the sampling time.

3. The imaging system as claimed in claim 2, wherein two-dimensional threads are allocated in the GPU, the thread size is 1024, and the block size is 4*fs*t / 1024.

4. The imaging system as claimed in claim 3, wherein the space size of the GPU video memory is 4*fs*t.

5. A photoacoustic remote sensing microscopy imaging method based on parallel processing acceleration, characterized in that: The imaging method comprises the following steps: Step 1: In the data stream mode, the acquisition card acquires PARS signals, trigger signals, and the position information of the fast axis and slow axis. The data of the four acquisition channels are in one data stream. The CPU creates a first buffer A and a second buffer B in the memory. The first buffer A is set as the receiving area, and the second buffer B is set as the processing area. The receiving area receives the data flowing out from the acquisition card, and if the processing area is not empty, the data in the area is copied to the GPU video memory. If the receiving area is full and the processing area has finished processing, the receiving area and the processing area are switched to continue receiving and processing data. Step 2: Allocate a two-dimensional thread in the GPU to perform effective data extraction on the sampled data; the extraction method is as follows: The index values ​​of the acquisition channel data at the same time are adjacent. The data of the four acquisition channels are divided with four data as a cycle. The data values ​​in the trigger signal data that are less than the threshold are set to zero. The threshold is set to 0.

3. The calculation formula is shown in (1): Calculate the rising edge position of the processed data, and the judgment formula is: Where f(n) is the original trigger signal data value, f′(n) is the processed trigger signal data value, n is the index value, and the index value that meets the above conditions is the position where the rising edge occurs; Step 3: Convert the valid data into the corresponding position data and PARS signal peak-to-peak data, and only calculate the valid data in the sampled data. The specific calculation method is as follows: After a single trigger, the average value of m valid data is calculated as the signal value at the time of single trigger to calculate the position data. The calculation formula is: Where X and Y represent the position information of the fast axis and the slow axis respectively, d is the focal length of the objective lens, k is the voltage to radian conversion ratio, sum x and sum y are the sum of m valid data of the fast axis and the slow axis, respectively, x and z y are the initial position voltages of the fast and slow axes, respectively; The peak-to-peak value is obtained by calculating the difference between the maximum and minimum values ​​of the subsequent m valid data of the PARS signal. The formula is: V PP =V max -V min (4) Where V PP is the peak-to-peak value, V max is the maximum value of the subsequent m valid data, V min It is the minimum value of the subsequent m valid data; Step 4: Convert the valid data into the corresponding position data and PARS signal peak-to-peak data. The extreme value seeking module obtains the extreme value of the position data based on the reduction algorithm of the alternating strategy, allocates two-dimensional threads in the GPU, and uses the GPU to calculate the position data according to the product of the thread size and the block size as the alternating factor to obtain an intermediate result array, and then copies the result array to the CPU; in the CPU, the final position data extreme value is calculated by serial processing; Step 5: Scale the X-axis and Y-axis position information to the final required total image pixel width. The calculation formula is: where x i and i are the X and Y coordinate values ​​of the i-th scanning point, respectively. max 、x min ,y max and min are the maximum and minimum values ​​of the X and Y coordinates obtained in the extreme value step; N x and N y are the number of grids that the image is divided into in the X and Y directions, i.e., the pixel width of the final image; n xi and n yi The grid coordinates corresponding to the i-th scanning point in the X and Y directions are respectively; according to the calculated coordinate information, a multi-atomic lock algorithm is used to apply for each grid's own atomic lock to synchronize the threads in the thread block. When multiple threads apply to access the same address, the atomic lock limits access to only one thread at a time, and the remaining threads wait for the atomic lock to be released; GPU statistics are used in the two-dimensional thread to obtain the number of scanning points contained in each grid, and the peak-to-peak value data in the valid data is used as the signal strength value of each scanning point. After each thread accesses, the result is visible to all threads; Step 6: Interpolate the grid. For grids containing multiple scanning points, calculate the average signal strength value according to the number of scanning points as the signal strength value of the grid; for grids that do not contain scanning points, use the linear interpolation method to fill them. The calculation formula is: Where i and j are the row and column coordinates of the grid, respectively. i+1,j 、V i,j+1 、V i-1,j 、V i,j-1 is the sum of the signal strength values ​​of each grid, N i+1,j 、N i,j+1 、N i-1,j 、N i,j-1 is the number of scan points of the respective grid, v i,j is the average signal value of the grid; the signal strength value data of all grid points calculated are used as grid data; Step 7: Use the grid data obtained in step 6 to find the extreme value using the reduction algorithm, and obtain the final image after grayscale normalization. The calculation formula is: Where i and j are the grid row and column coordinates, I i,j is the pixel value after grayscale normalization, v i,j is the grid signal value, v max 、v min are the maximum and minimum signal strength values ​​of all grids.

6. The imaging method according to claim 5, wherein the space sizes of the first buffer A and the second buffer B are both 4*fs*t, wherein fs represents the sampling frequency of the acquisition card, and t represents the sampling time. 7 . The imaging method according to claim 6 , wherein the two-dimensional threads are allocated in the GPU, the thread size is 1024, and the block size is 4*fs*t / 1024.

8. The imaging method according to claim 7, wherein the space size of the GPU video memory is 4*fs*t.

Images