Implementation of 2-Dimensional Medical Image Reconstruction Algorithm for Photoacoustic Imaging

In this paper, we review and discuss our proposed image reconstruction algorithm for two dimensional stereo photoacoustic medical imaging system. The proposed algorithms aim to evaluate the existence possibility of acoustic source within a search area by using the geometric information, which is the distance between each sensor element of ultrasound probe and corresponding searching point. Besides, we implemented the hardware and software aspects of a photo-acoustic imaging system that described in hardware implementation section of this paper. At the end of this paper, some simulation results are presented which utilized the proposed algorithm while compared the image quality of its and conventional k-wave algorithm in which FFT should be used. The simulation results proved the effectiveness of the proposed reconstruction algorithm. Keywords— Photo-acoustic; k-wave; medical imaging; GPU; NON-FFT; reconstruction algorithm


I. INTRODUCTION
Photo-acoustic tomography is an emerging biomedical imaging modality. It is useful for visualizing a structure embedded in soft tissue, which absorbs light [2]. It based on the generation of ultrasound waves within a structure embedded in soft tissue due to the absorption of short pulses of laser light. When a laser pulse illuminated upon a structure embedded in soft tissues, the resulting heat related expansion of the structure in soft tissues produces a local increase in pressure and this increase in pressure propagates outward as an ultrasound waves. Depending on the optical wavelength of the input laser light, photons absorbed by endogenous chromosphere such as water and hemoglobin and heat transfer due to optical absorption cause a thermoelastic expansion of a structure in embedded soft tissue [3]. This thermoelastic expansion makes increase in acoustic pressure and this makes generating acoustic wave on the surface of structure. The acoustic waves generated are usually broadband with energy at frequencies from hundreds of kilohertz to many tens of megahertz. By analyzing the acoustic waves, we can find the information of distribution of laser beam absorption within a region of soft tissue. This distribution information, which are associated with cellular level of physiological status such as normal or tumor, can transformed into contrast in imaging. Molecular imaging is also possible by using targeted exogenous contrast agent [4]. Large absorption difference between normal and abnormal tissue makes it promising to use laser beam as a thermo-acoustic tomography to detect cancers. The acoustic waves arriving at the tissue surface are detected using one dimensional or two-dimensional array of ultrasound detectors. Using the recorded ultrasound time signals, an image of the initial pressure distribution can be reconstructed. The reconstruction of photo-acoustic image requires a mapping form the time varying ultrasound signal recorded on the tissue surface to the three-dimensional initial pressure distribution in side of tissue. A large number of reconstruction algorithms are currently available. Almost all algorithm based on inversion of Radon transform. Many results of Radon inversion algorithms are developed and applied to the research field of radiation oncology [5]- [8]. The most common algorithm used in X-ray CT is calling filtered back-projection [9]- [15], [19]- [21]. While mathematically explained perfectly, computational load is too heavy to implement by using software and this fact makes restriction in their practical usage. Another algorithm based on the harmonic decomposition of the imaging problem have been proposed. For a planar detector, this algorithm can implemented by using FFT and inverse FFT. Time-reversal algorithm is another method in image reconstruction [17], [18]. In this algorithm, image reconstruction achieved by running a numerical model of the forward problem backwards in time. Nevertheless, if there exists nonlinear ~ 7 ~ characteristics in the propagation of acoustic wave, this algorithm degrade the quality of the reconstruction image. Recently, this time-consuming part can be handled and solved by using multithreaded parallel programing due to the help of technology advance in graphic processing unit (GPU). GPU and FPGA (Field Programmable Gate Array) hardware platforms are hot research issues in the field of medical imaging since parallel processing speed is ten or more times faster than multi-core CPU hardware platforms [16], [22], [23]. In this paper, we review and discuss again our proposed image reconstruction algorithm for two-dimensional stereo photo-acoustic medical imaging system. The proposed algorithms aim to evaluate the existence possibility of acoustic source within a search area by using the geometric information which is the distance between each sensor element of ultrasound probe and corresponding searching point. Besides, we implemented the hardware and software aspects of an photo-acoustic imaging system that is described in hardware implementation section of this paper. At the end of this paper, some simulation results are presented which utilized the proposed algorithm while compared the image quality of its and conventional k-wave algorithm in which FFT should be used. The simulation results proved the effectiveness of the proposed reconstruction algorithm. Both FPGA and GPU were used in order to develop and implement a simple and fast photo-acoustic medical imaging system based on the proposed algorithm shown in Fig.1.

II. NEW RECONSTRUCTION ALGORITHM
In our previous paper, we proposed a new image reconstruction algorithm based on the geometric information of the search area and acoustic sensor under the assumption that the speed of acoustic wave is constant. Also, the paper intelligibly explained proposed algorithm with 1D and 2D sensor array cases [1]. In this section, we briefly describe the equation of the proposed algorithm with the calculation for arc length algorithm completely described. First, let us define the geometric parameters as shown in Fig.2.
where i 2 [1;Ny]; j 2 [1;Nx]; k 2 [1;Ns], where p(k; t) is the acoustic signal of k-th sensor element, T is the sampling time of the acoustic sensor signal, t(i; j; k) denotes the corresponding time delay for d(i; j; k), l(i; j; k) is an arc length whose radius is d(i; j; k) in the searching area, d(i; j; k) is the distance between the searching area corresponding to the position indexed by p(i; j) and k-th element of acoustic detector, m0(i; j) is total measure of the strength for the acoustic source at location indexed by p(i; j), and α is any positive constant chosen considering the attenuation factor of the acoustic wave in a tissue. Arc length calculation algorithm is described a little bit in our previous paper and it did not use any simulation. Next subsection completely describes how to implement the algorithm and gives some results.

A. Research Goal
In this subsection, we explain how to calculate a part of arc shown in Fig.3 which is needed in implementation of the proposed algorithm with arc length calculation. The general idea is to find the hitting points and then calculate an arc angle and length of the arc. At first, find the hitting point of circle shown in Fig.3 and the four extended straight lines of the rectangle which marks the searching area of acoustic sources. In order to find the hitting points in the search area, we have been considering four conditions such as TOP, LEFT, RIGHT, and BOTTOM side of the searching area. In the Fig.3, the line sensor is center of a circle. Therefore, a maximum of one point can be found for the left and right sides while a maximum of two points can be found for the top and bottom sides. The following equations describe how to find a hitting point in all side of the searching area. First, find a hitting point in the top side. If one or two hitting points are found, Y-axis coordinator of the hitting point is equal to yt. If X-axis coordinator of the hitting point is between xl and xr, the hitting point is found in top side of the search area. The hitting points of top side can be expressed by using the following equation of circle.
If (( − ) 2 + 2 − 2 )is negative, there exists no hitting point. Second, find a hitting point in the left side. If a hitting point is found, X-axis coordinator of the hitting point is equal to . If Y-axis coordinator of the hitting point is between and , the hitting point is found in left side of the search area. The If (( − ) 2 + 2 − ( − ) 2 ) is negative, there exists no hitting point. Third, find a hitting point in the right side. if a hitting point is found, X-axis coordinator of the hitting point is equal to xr. If Y-axis coordinator of the hitting point is between the yt and yb, the hitting point is found in right side of the search area. The hitting point of right side can be expressed by using the following equation of circle.
In this case, we can find four different hitting points (one, two, three and four). Fig.4 shows that how many hitting points can be detected in this simulation. The arc length is zero if the number of hitting points is one. If there are two hitting points, then calculate the angle between these two points and obtain the arc length by using angle and radius of the arc. If there three hitting points, one of the three hitting points is in left or right side and another one is in bottom side of the search area. The last hitting point is equal to one of the bottom corner point in search area. In other words, an arc that is outside of the search area, between the corner point and another bottom point. Therefore, we do not need to calculate the arc and remove the hitting point in the bottom corner. Now, there are only two hitting points and we can then calculate an arc using the previous algorithm for two hitting points. If four hitting points are found, two of them are in left and right side while another two are in bottom side of the search area. Moreover, two different arcs must be calculated, one of the arc is calculated by using the hitting point in left side and the hitting point of left in bottom side. Another arc is calculated by using the hitting point in right side and the hitting point of right in bottom side.

III. HARDWARE IMPLEMENTATION
In this section, the hardware solution of simple medical ultrasound system and how to apply the proposed the algorithm is introduced. In order to develop and implement photo-acoustic medical imaging system which is exactly appropriate for the proposed algorithm, both FPGA and GPU are used. The signal receiving part of this system is initially tested using the 8channels AFE5808A-EVM and TSW1400EVM which is shown in Fig.5. Even though the evaluation boards are working fine, the VCA (Voltage Controlled Attenuator, included in the AFE5808A) cannot be fast enough handled via a digital interface. Therefore, AFE5816-EVM has 16 channels and its attenuator can be controlled fast enough by digital signal while the output of the ADC (Analog to Digital Converter) is a lowvoltage differential signaling (LVDS) [25], used instead of AFE5808AEVM. Fig.6 shows environment that implementation of the development system with 16 channels AFE module, L14-5/38 ultrasound probe, and TSW1400EVM captured board. Fig.7 and Fig.8 show the timing diagram of raw signals from AFE5808A and AFE5816 devices by using TSW1400EVM while see that controlled time of TGC. From these hardware experiments, AFE5816 device attenuator can be handled swiftly via digital interface. In other word, we can consider TGC strategy into the hardware section of AFE (Analog-Front-End) instead of software section. The next section describes and shows that why TGC method is critical in medical image reconstruction algorithms, especially the proposed algorithm.

IV. SIMULATION RESULTS
Two cases were simulated in order to check and prove the proposed algorithm is better than conventional k-wave algorithm. These simulations, one of the k-wave examples that uses a line detector is chosen for the following numerical illustration, while the sonar speed is defined as 1500m=sec in these simulations. In first simulation, the five acoustic sources are respectively located at considered whose radius and magnitude are set to 1mm and 5, respectively. In second simulation, three acoustic sources are located anywhere in the region, 0 ≤ ; ; ≤ = 100: 0 . Three acoustic sources were considered whose radius and magnitude are set to 5mm and 1, respectively. Fig.7 (a) and Fig.9 (a) show the locations and sizes of the original acoustic sources while Fig.7 (b) and Fig.9 (b) are illustrated raw signals of some elements of linear detector. In the first simulation, five acoustic sources are used even though a sensor can be detect five or less peak local signals. In other word, a signal of two or more acoustic sources can be added each other. For instance, in Fig.7 (b), the second raw signal of 25th element of acoustic detector has four peak magnitudes. Fig.8 (a) and Fig.10 (a) gives the image reconstructions obtained by using conventional method in which FFT and IFFT transform is used, while Fig.8 (b) and Fig.10 (b) show the simulation results obtain by using the proposed algorithm. The images qualities such as shape and magnitude of the proposed algorithm is better than k-wave conventional algorithm. From these two examples, the image quality such as shape got worse and the conventional algorithm is barely able detect a source that is located far away from linear detector. Because the proposed algorithm considered TGC method, its is clearly detected acoustic sources which is located in anywhere in the searching area.

A. Conclusions
In this paper, the implementation of hardware and software our proposed image reconstruction algorithm for two dimensional stereo photo-acoustic medical imaging system. The hardware setup for the photo-acoustic tomography system and some experimental results of the system are also given. The hardware implementation section describes how the proposed algorithm and hardware setup are combined with TGC method based on AFE device. Simulation results using Matlab k-wave toolkit illustrate the effectiveness of the proposed algorithm. In addition, the quality of reconstructed image is better by using the proposed algorithm compared to the convolution of k-wave algorithm.

B. Future Works
We are working on a project with Samsung Medison company for developing medical photo-acoustic instrument for cardiovascular and plaque of carotid artery imaging. Nowadays, many researchers are trying to implement image reconstruction software by using both FPGA and GPU hardware technology. In our laboratory, an environment through a new embedded system is developed and implemented. It is consists of NVIDIA Tegra K1 embedded board which is based on GPU technology and Kintex-7 FPGA KC-705 Evaluation board for medical photo-acoustic system to reduce the calculation time. Furthermore, we have been developing reconstruction algorithm in Matlab with k-wave toolbox and on NVIDIA Quadro K4200