1. Introduction

Prostate cancer (PCa) is one of the most common malignant tumors among men worldwide, with a high mortality rate [1]. According to [2], 11.6% of men will develop prostate cancer in their lifetime, with approximately a 20% death rate in the United States. Radical prostatectomy is a popular surgical approach to treat PCa by removing the entire prostate gland since 1905 [3,4]. In clinical practice, the traditional open radical prostatectomy (ORP) has almost been replaced by laparoscopic radical prostatectomy (RLP) [5]. As a minimally invasive surgical procedure for PCa, RLP significantly reduces blood loss, hospitalization duration, and postoperative complications [6]. However, the long learning curve associated with laparoscopic procedures limits the application of RLP [7]. Robot-assisted laparoscopic prostatectomy (RALP) has been demonstrated [5] to shorten this learning curve by leveraging the wristed instruments and the 3-D endoscopic camera of the telerobotic surgical system, usually the da Vinci surgical system, to achieve intuitive operation [8]. However, the endoscopic camera cannot localize the prostate lesions nor visualize the sub-surface anatomy of the prostate gland. Therefore, a complementary medical imaging modality is necessary to facilitate RALP.

Transrectal ultrasound (TRUS) facilitates the intraoperative identification of vital anatomic structures of the prostate gland and adjacent tissues, e.g., peripheral nerve and neurovascular bundles (NVBs), in a non-invasive and real-time manner [9,35]. TRUS provides precise ultrasound images of tissue pointed by the surgeon with the instrument. Manipulating the TRUS is challenging during RALP, because the da Vinci surgical system is positioned between the patient's legs, thus restricting access to the TRUS [10]. Furthermore, the surgeon seated at the console must rely on the presence of a dedicated ultrasound assistant to operate the TRUS probe. To alleviate the challenges associated with manual manipulation of the TRUS, a robotic transrectal ultrasound system for prostatectomy was presented and has been adopted in RALP studies [10–12]. The da Vinci and TRUS coordinate systems are registered by touching the tissue with the instrument tip at three or more different locations [30]. However, the search process for each location is manual, time-consuming, and requires an additional skilled operator in the operating room [11]. Besides, surgical tool tip detection in ultrasound images is challenging, especially in the elevational direction, because it is dependent on the shape of the surgical instrument and the way it is pressed on the tissue [19,30].

A very accurate alternative registration method between the TRUS and the da Vinci system involves the use of photoacoustic markers (PMs), which are generated by the photoacoustic (PA) effect, i.e., the formation of acoustic waves when a pulsed laser source is absorbed by optical materials. Robot assisted photoacoustic imaging using the da Vinci surgical system was investigated in [13–16]. In [17], authors demonstrated that photoacoustic markers generated by a pulsed laser could be reconstructed as images representing the local optical absorption of tissue constituents. Cheng et al. [18]’s research proved the feasibility and advantages of PM in surgical guidance. Instead of physically touching the tissue, a PM-based approach using focused laser illumination as the PM was introduced in [19] for RALP. In this method, the endoscopic stereo camera frame was registered with respect to the TRUS frame using Horn’s method [33]. While accurate, this registration is still time-consuming (as much as 7 min), because the TRUS is manually rotated to align its imaging plane to the laser spots. Driven by this limitation, an automatic image-based search for the photoacoustic markers is proposed in this paper.

The directivity pattern of clinical transducer was investigated in [20]. The photoacoustic image intensity is proportional to the distance between the PM and the TRUS imaging plane. When the PM is coplanar with the imaging plane of the TRUS transducer, the detected PA image intensity of the PM is maximum [21]. Therefore, we have a PA image intensity function, which, for a single PM, should be unimodal and symmetric with respect to TRUS rotation [20]. Based on this, the key component of our work is to conduct a one-dimensional search in the presence of PMs to find the TRUS rotation angle corresponding to the maximum PA intensity.

The numerical optimization of unimodal functions in one-dimensional (1-D) search space is a common topic in many engineering fields [22–24]. Introducing an efficient algorithm to replace the exhaustive search will significantly shorten the duration of the registration. The Golden Section Search (GSS) algorithm [25] is a classic approach to solve this maximization problem. Similarly, Ternary Search (TS) algorithm [27] is a widely used divide-and-conquer algorithm for searching the maximum of a unimodal function [25,26].

In this paper, we propose a novel automatic search module to localize the PMs using the robotized transrectal ultrasound actuator in robot-assisted radical laparoscopic prostatectomy. The key advance presented in this paper is an efficient manner to replace the time-consuming manual search required to register the TRUS to the robot coordinate systems. We simulate and experimentally evaluate several search algorithms. The proposed framework outperforms the manual search in terms of time and does not require additional staff in the operating room.

2. Method

2.1 Problem definition

In our RALP system depicted in Fig. 1. (a), a laser source generates a PM on the surface of the tissue. The PM is in the plane of the TRUS when the TRUS angle is ${\theta _p}$. The PM can be detected by a TRUS with a linear array of the Ultrasonix biplane TRUS brachytherapy probe, which has a width of 55 mm, consists of 128 elements, and operates within a 5-9 MHz bandwidth. The position of TRUS imaging plane is denoted as ${\theta _i}$. After delay-and-sum (DAS) beamforming, the PA image is used to calculate the PM intensity. The detected intensity of the PM achieves a maximum when the TRUS imaging plane is coplanar with the PM, namely the TRUS position ${\theta _i}$ equals ${\theta _p}$, and decreases when the ${\theta _i}\; $ moves away from ${\theta _p}$. As described in Fig. 1. (b), there is a PA signal intensity function related to the rotation angle of the TRUS. Theoretically, for a single PM, this intensity function is unimodal and symmetric [20]. Our goal is to search for its peak. We can drive the motorized TRUS to rotate within -35° to 35°, providing enough coverage for the prostate imaging, and synchronously calculate PA signal intensity at each location using the US + PA imaging system. In this way, we will investigate efficient algorithms to solve this problem compared to manual scanning.

 

Fig. 1. (a) The laser illuminates a small spot (PM). The motorized TRUS is located at ${\theta _i}$. The goal is to find ${\theta _p}$ of the PA marker. (b) The TRUS also receives out-of-plane signals. Therefore, PA signals will be recorded at different angles.

Download Full Size | PDF

2.2 Weighted average

The laser illuminates a small spot as shown in Fig. 1. (a). The motorized TRUS rotates from -35° to 35° in equal increments of ${\theta _d}$ at ${\theta _i}$s where $i = 1,\; \ldots ,N$ for a full scan. The goal is to find ${\theta _p}$ in the TRUS coordinate system. The TRUS also receives out-of-plane signals. The maximum PA signal intensity, $S({{\theta_i}} )$, is recorded for each ${\theta _i}$, described in Fig. 1. (b). As the directivity pattern of the TRUS transducer with respect to the incoming waves is symmetric, a simple weighted average can be used to localize the PA marker, ${\theta _p}$:

Note that we need to set a threshold to filter out the samples with very low or no PA signal intensity. This helps to reduce noise and improve localization accuracy.

2.3 Golden section search and ternary search

The Golden Section Search (GSS) is an efficient algorithm for finding the extremum (maximum or minimum) of a unimodal function. The GSS has been widely used in many applications because of its rapid convergence, without requiring that the objective function be continuous. In the GSS algorithm, we compare a pair of observation at two interior points ${\theta _{m,1}}$, ${\theta _{m,2}}$ ${\in} R\; $ satisfying ${\theta _l} < {\theta _{m,1}} < {\theta _{m,2}} < {\theta _r}$ for each iteration. The two interior points partition the search interval in terms of the golden section ratio r. Here, we set the initial value of ${\theta _l}$ and ${\theta _r}$ as -35° and 35°. They are updated at each iteration and eventually convergence towards ${\theta _p}$. The interior points ${\theta _{m,\; 1}}$ and ${\theta _{m,2}}$ are computed as follows:

The strategy of GSS can be stated as follows: if $S({{\theta_{m,1}}} )\ge S({{\theta_{m,2}}} )$, the maxima must be in the interval $[{{\theta_l},{\theta_{m,\; 2}}} ]$. Likewise, if $S({{\theta_{m,1}}} )< S({{\theta_{m,2}}} )$, the maxima must be in the interval $[{{\theta_{m,1}},{\theta_r}} ]$. Subsequently, we compute two new interior points in the interval mentioned above for the next iteration. At every iteration, we calculate the interior points in this manner and keep narrowing the search interval until meeting the termination condition. Theoretically, the iteration would terminate when ${\theta _l} \ge {\theta _r}\; $ is met. In practice, however, we can set a termination condition to boost the speed of convergence. Once the search interval is less than this termination condition, we stop the iteration and take ${\theta _p} = ({\theta _l} + {\theta _r})/2$ as the output location. If we set r = 1/3 instead of the golden section ratio of 0.618, the algorithm turns into the Ternary Search (TS). Except for the scalar r used for interior points calculation, the other parts of the TS are the same as the GSS.

2.4 Evaluation method

The proposed automatic search module was validated by both simulation and experiment (including offline analysis and online analysis). The performance of various search algorithms was evaluated in terms of localization accuracy and time efficiency. For the evaluation of localization accuracy, we measured the localization error which was defined as the absolute value of the difference between the reference location and optimal output location. To assess time efficiency, we adopt the number of data acquisitions as the metric due to most of the time in the proposed search procedure being used to acquire data.

2.5 Simulation validation

For the simulation, the experimental directivity pattern and PA model developed in [20] was applied to the point source PA signals generated by MATLAB k-Wave toolbox. The simulated receiving transducer models a standard linear array of Ultrasonix Biplane TRUS transducer, which has 128 sensing elements in 55 mm width. Its center frequency is 7.5 MHz. The grid size in k-Wave simulation is 0.1 mm ${\times} $ 0.1 mm ${\times} $ 0.1 mm. The intensity of simulated PA signals was normalized to a value between 0-1. The reference angles of PA markers were located at 1000 random locations within -35$^\circ $ to 35$^\circ $, and PA signal was recorded at the desired positions based on the proposed search algorithms. Here, we compared the performance of all three search algorithms with different numbers of data acquisitions. In the GSS and TS algorithms, the angle of the TRUS was determined by the policy introduced in 2.3 and the iteration/acquisition could be continued until the desired accuracy was met. The number of data acquisitions was a multiple of two in these two algorithms, because each iteration needed to acquire PA signals twice. We tested the performance of these two algorithms with 2 to 28 acquisitions. With the WA algorithm, laser spots were located at a random position and PA signals were detected at several ${\theta _i}$s spaced at ${\theta _d}$ apart, with

So, as the number of acquisitions, N, increases, ${\theta _d}$ decreases. Note that the number of acquisitions could not be too small because there may be no PA signals generated, which is why we started with a minimum number of 8 acquisitions. The number of acquisitions was a multiple of two to keep it consistent with the GSS and TS algorithms. Thus, for the WA algorithm, we studied the performance when acquiring data 8-28 times.

To investigate the effect of thresholding the PA signals in the WA algorithm, we ran the WA algorithm simulation with various thresholds between 0 and 0.35 with an interval of 0.05.

2.6 Experiment validation

In addition to the simulation validation, we conducted experimental studies to illustrate the advantages of our proposed method. In our experiment, a pulsed laser diode illuminated random locations on the tissue surface. To validate our approach, we conducted two types of analysis: offline and online. For the offline analysis, we conducted experiments to collect data by scanning the entire rotation interval with a very small step size (∼0.5$^\circ $), enabling us to have a bank of PA data for our investigation. We simulated the online search based on this collected data. The numerical precision of the motor rotation angle was made consistent with the online analysis by interpolating the original offline data. For the online analysis, the algorithms were investigated in a real-time implementation where each of the search algorithms were integrated into the robotic and imaging systems. As a result, PA data were acquired at the locations that the algorithm asked for, and the output was compared with the reference acquired from a manual search.

A. System architecture

The overall system depicted in Fig. 2 comprises of several modules including the main host computer, pulsed laser system, TRUS actuator control module, and PA + US imaging system. The imaging system was described in [19]. Here, we emphasize the modules related to the automatic implementation of the search algorithms. In the pulse laser system, the 785 nm pulse laser diode (QSP-785-4, QPhotonics LLC, Ann Arbor, MI, USA) is excited by a laser driver (LDP-V, PicoLAS GmbH, Germany) to generate a 1µJ/pulse laser beam with a 300 µm diameter. The PA + US imaging module synchronously acquires PA signals. The generated PA signals can be detected by a TRUS transducer (BPL 9-5/55, BK Medical, Peabody, MA, USA) and a SonixDAQ (Ultrasonix Medical Corp., Canada). The actuator control module can drive the TRUS transducer to rotate through an actuator (rotation stage, PRM1Z8, Thorlabs, Newton, NJ, USA) rigidly attached to it. Real-time communication between each component in the system is based on the Robot Operating System (ROS).

 

Fig. 2. Schematic diagram of the overall system composition. It also describes the TRUS coordinate system (TR). The y-axis is parallel with the axis of rotation of the TRUS.

Download Full Size | PDF

B. Experimental setup

The experimental set-up is shown in Fig. 3. (a). An ex vivo bovine liver tissue was supported by two pieces of ultrasound gel pad (Parker Laboratories, Fairfield, NJ, USA) with a total thickness of 5 cm to mimic the prostate gland. A laser diode held by a clip was placed close to the surface of the ex vivo tissue. The laser's location and its surrounding region (5mm${\times} $5 mm area) were stained with India ink (Dr. Ph. Martin, Oceanside, CA, USA) to increase laser absorption. A TRUS transducer motorized by an actuator was placed beneath the ex vivo tissue and ultrasound pads to detect the PM. Here, we tested the proposed approach on a homogenous black plastisol phantom, and three types of ex vivo tissue, bovine liver, muscle, and chicken breast, to demonstrate its feasibility. In actual surgical scenarios, exogenous optical contrast agents, such as indocyanine green (ICG) and voltage-sensitive dye (VSD), are used for fluorescence endoscopic imaging [36,37]. Thus, we utilized ICG dye (TCI, Tokyo, Japan) to stain ex vivo bovine muscle tissue in one experiment, aligning with the protocols for prostate nerve sensing [28], to evaluate the effectiveness of the proposed methodology in a practical surgical context. As illustrated in Fig. 3. (b), the PA signals arranged on a 3${\times} $3 grid were collected on each type of material. The distribution of PM occupied most of the space in the TRUS field of view. In addition to the PA sensitive range of the transducer, the TRUS actuator can rotate to any position within -35° to 35°, with a minimum incremental motion of 25 arcsec (0.0001 degree). Thus, there is no blind area for detecting photoacoustic signals. Note that the initial position of TRUS was always at -35°.

 

Fig. 3. (a) The offline analysis experimental set-up with bovine liver ex vivo tissue stained with India ink. (b) The conceptual illustration of the laser spot location on ex vivo chicken breast tissue from the top view of the x-y plane. The red dots indicate the laser spot, and a 3 × 3 grid roughly spaced at 20 mm in the x-axis and y-axis directions, respectively, of the x-y plane of the TRUS coordinate system. For other ex vivo tissues and plastisol phantom, the laser spot location is arranged in the same way.

Download Full Size | PDF

C. Online search implementation

We integrated the automatic search module into the overall TRUS-da Vinci robotic system using ROS and used a LabVIEW script [19] to control the rotation stage. To integrate this TRUS controller into our ROS communication system, we used the ROS for LabVIEW toolkit software [29] developed by the Tufts University Center for Engineering Education and Outreach. We created a node for the TRUS controller, as well as a publisher and a subscriber to receive the angle that the TRUS should rotate to and return the actual rotated angle. For online analysis, there was no need to exhaustively scan the entire rotation range with a small step size to collect data as in our offline analysis. The TRUS only rotated to specified locations that the algorithm asked for and acquired PA signals from the data acquisition system at these locations. Here, we tested the proposed online automatic search module on a black plastisol phantom. Like in the offline analysis case, we tested each search algorithm on the 3${\times} $3 grid of locations and the results were compared with manually found reference locations.

D. Tool Tip Detection - Sensitive Range

When using the da Vinci instrument tool tip for registration, a key stage in this procedure is the da Vinci surgical instrument tip detection [30]. Similar to the PM intensity, the tip intensity is the highest when the US imaging plane is coplanar with the surgical instrument tip. It decreases and eventually disappears as the US imaging plane deviates from the coplanar position. According to the previous work [30], for registration, we only need to make sure the surgical instrument tip is within the range (termed sensitive range) that can be detected by the US, not necessarily to be exactly coplanar with the US imaging plane. The role of tool tip detection is similar to the PM search for registration. Thus, we would like to evaluate the search error by comparing it to the sensitive range of tool tip detection.

We experimentally investigated the sensitive range when detecting the da Vinci surgical instrument tip. As shown in Fig. 4., the da Vinci instrument tip held on a robotic arm was pressed against the surface of a piece of ex vivo chicken breast. We rotated the TRUS clockwise with a step size of 1$^\circ $ starting with the reference location and recorded the corresponding B-mode image at each location. The range from the reference location to the location where the tip disappeared (shown in Fig. 10. (b) and (c)) was half of the sensitive range. Assuming that the function of detected tip intensity versus TRUS rotation location is symmetric, the other half sensitive range was inferred by flipping along the reference location. Here, we roughly calculated the entire sensitive range by adding these two half sensitive ranges. Three locations of the da Vinci surgical instrument tip spaced at about 2 cm in x-axis were tested in this experiment.

 

Fig. 4. The experimental setup for investigating the sensitive range of da Vinci surgical instrument tip detection.

Download Full Size | PDF

For the B-mode image processing, we first segmented the da Vinci instrument tip by leveraging GrabCut [31], following to find the centroid of this connected instrument tip area, and finally taking the pixel value of this centroid as the instrument tip intensity. Let ${\theta _{high}}$ denote the location with the highest tip intensity, and ${\theta _{low}}$ denote the location with the lowest tip intensity that can be detected. Thus, half of the sensitive range can be defined as the absolute value of these two locations, namely $|{{\theta_{high}} - {\theta_{low}}} |$.

3. Results

3.1 Simulation results

For the TS and GSS algorithm, we tested a set number of acquisitions from 2 to 28 in increments of 2. Similarly, we tested a set number of acquisitions from 8 to 28 for the WA algorithm to ensure that at least one valid PA signal can be acquired at the proposed equal intervals. We investigated the effect of the threshold for the WA algorithm by running the simulation under various thresholds. As depicted in Fig. 5., the average final error dropped as the threshold rose in the range of 0-0.25. But after exceeding this range, the error increased since meaningful signals were also removed by thresholding. We chose a threshold of 0.1 in the WA algorithm, because it provided a low error across the entire range of the number of acquisitions. Note that even without a threshold, we still obtained a small average error of 1.02°. The average error after 1000 trials changed with the number of acquisitions as shown in Fig. 6. The error of TS and GSS algorithms rapidly dropped at the beginning and converged to a small value, while the error of WA did not change much with the number of acquisitions. When conducting 28 acquisitions, the average errors of the TS and GSS algorithm were 0.06$^\circ $ and 0.01$^\circ $, respectively. The GSS and TS algorithm gradually approach the target with the number of iterations. The simulation results were consistent with this property. In the WA algorithm, the average error when acquiring different numbers of PA signals was 0.47$^\circ{\pm} $0.11$^\circ $.

 

Fig. 5. The average error of 1000 experiments using the WA algorithm with different thresholds in the simulation study.

Download Full Size | PDF

 

Fig. 6. The average error of 1000 experiments using various algorithms in the simulation study. Red dots, blue dots, and green dots represent the error of TS, GSS, and WA algorithm.

Download Full Size | PDF

3.2 Offline analysis results

In our offline analysis, we interpolated each dataset to achieve the same resolution in the rotation angle as the simulated data. Figure 7 depicts the average error versus the number of acquisitions on various materials resulting from the offline analysis. Table 1 summarizes the error of each algorithm when acquiring the PA signal 16 times. Note that we only obtained one set of data on bovine muscle stained by ICG, so the results of this test did not have a standard deviation.

 

Fig. 7. The average error of offline analysis versus the number of acquisitions. The sub-figures (a) - (e) indicate the offline analysis results on bovine liver, muscle, chicken breast, phantom and bovine muscle stained with ICG, respectively.

Download Full Size | PDF

Table 1. Average Error of Various Algorithms in Offline Analysis Terminating After 16 Times of Acquisitions

View Table | View all tables in this article

3.3 Online analysis results

In online analysis, each algorithm described earlier was integrated into the TRUS-da Vinci robotic system. When conducting the online analysis, the TRUS automatically rotated to the position specified by the algorithm, and finally stopped at ${\theta _p}$ calculated by this algorithm. Figure 8 illustrates the comparison of the average errors of TS and GSS algorithms with varying numbers of acquisitions. The average errors of these two algorithms decrease with an increase in the number of acquisitions. For the TS and GSS algorithms, the specific location of the TRUS at each number of acquisitions is illustrated in Fig. 9. The GSS and TS algorithm achieved minimum error, 0.48°${\pm} $0.30$^\circ $ and 0.29°${\pm} $0.30$^\circ $, when acquiring 28 times data. When adopting WA, we took 5$^\circ $ as the step size, that is, acquiring 15 times PA signals to cover the entire space. However, in the online analysis, we can terminate the search as long as the entire peak of the intensity is recorded and passed. We can often complete the WA search by collecting data fewer than 15 times. The average error of WA was 0.53°${\pm} $0.30^o with 8.7 data acquisitions. Table 2 is the comparison of the average error of various algorithms when taking a similar number of acquisitions.

 

Fig. 8. Average error of TS and GSS with a varying number of acquisitions in online analysis. The blue bar and orange bar indicate the TS and GSS algorithms respectively.

Download Full Size | PDF

 

Fig. 9. TRUS location during the TS and GSS algorithms. The blue dots, red dots, and orange dashed line refer to the locations in TS and GSS as well as the reference location ${\theta _{gt}}$. The position 1 to 9 were located on a 3${\times} $3 grid as described in Fig. 3. (b).

Download Full Size | PDF

Table 2. Average Error of Various Algorithms in Online Analysis

View Table | View all tables in this article

3.4 Sensitive range results

We pressed the da Vinci surgical instrument tip against the ex vivo chicken breast tissue surface at three locations. For each tip location, a set of B-mode images were recorded with the TRUS rotated in steps of 1$^\circ $ (Fig. 10. (a) and (d) show one of the B-mode images sets and the corresponding tip intensity curve). We took the position corresponding to the strongest and weakest tip intensity that can be detected as half of the sensitive range. In our experiment, half of the sensitive range of the three locations were 5$^\circ $, 4$^\circ $, and 4$^\circ $ respectively, and the average value was 4.33$^\circ{\pm} $0.58$^\circ $.

 

Fig. 10. (a) A set of B-mode images with varying tip intensity. The da Vinci surgical instrument tip is within the red dashed circle. The reference location is 0$^\circ $. (b) The B-mode image at the reference location. It can be partitioned into 3 layers: from top to bottom are air, ex vivo tissue and ultrasound pad. (c) The B-mode image at a location transcending the sensitive range. (d) The tip intensity versus TRUS location. As the TRUS location moves away from the reference location 0$^\circ $, the corresponding tip intensity keeps reducing.

Download Full Size | PDF

4. Discussion

Sensitive range of tool tip detection is a practical way to evaluate the search errors of the proposed methods. In terms of the experimental results described in 3.4, the sensitive range for the da Vinci surgical instrument tip detection is approximately 8.66$^\circ $. In fact, we can obtain a much lower error (0.53$^\circ{\pm} $0.30$^\circ $) using the WA algorithm in our study just acquiring PA signals eight times. Consequently, the accuracy of all three algorithms should be considered precise and sufficient compared to the sensitive range of tip detection. Regarding the three algorithms adopted in this paper, TS and GSS output the most precise PM location on ex vivo tissues when acquiring data 28 times.

The proposed search algorithms have the potential to be used as a part of the control platform for the surgical tool tracking and TRUS servoing tasks. In online analysis, and with only acquiring two PA signals, the average error was 4.20$^\circ $ and 5.89$^\circ $ when adopting the TS and GSS algorithm starting at the far left (-35°) and right (+35°) of the boundary which is reasonable when compared to the required 8.66$^\circ $. However, this accuracy will significantly improve with more data acquisitions as illustrated in this paper. Besides, the TS and GSS tend to converge toward and settle at the target location ${\theta _p}$. This property makes these two algorithms potentially applicable in tool tracking procedures and will be further investigated in our future studies.

The time complexity of the WA algorithm is $O(n )$, while the time complexity of the GSS and TS algorithms is $O({{{\log }_3}n} )$. So, the GSS and TS algorithms take less time theoretically. In practice, this computational time is negligible relative to the time spent on reading data from the SonixDAQ and rotating the TRUS. Therefore, we investigated the number of data acquisitions and the travel angle of the TRUS actuator as indicators of the required time. In this work, we utilized the API of SonixDAQ to read raw PA signals from the device and read these raw data files to a PC for subsequent beamforming and reconstruction. This data transfer typically takes approximately 10 seconds. Driven by this limitation, our group is developing a real-time PA imaging based on SPARE (Synthetic-aperture Based Photoacoustic Re-beamforming [32]). SPARE is a new PA reconstruction approach, which uses post-beamformed radio frequency (RF) data from an ultrasound machine as the input data rather than the raw channel data in the DAQ. Since no data transfer time from DAQ to PC is required, real-time PA imaging based on SPARE is very possible. Real-time PA imaging can essentially benefit the time efficiency of the PM-based RAPL system. Besides, improving PA imaging speed [38] is another solution. In that case, the only time we need to consider is the time of TRUS rotation. When acquiring data for 16 times, the total rotation travel in TS and GSS are about 36$^\circ $ and 24$^\circ $, respectively. Thus, the estimated time of rotation is within 1-2 seconds, which is much less than the 2-3 minutes of manual search reported in [30] and close to a real-time search approach.

Regarding the selection of tissues for the ex vivo study, obtaining Institutional Review Board approval for the use of ex vivo human prostate tissue is challenging and requires preliminary results. Therefore, in this study, we opted for three readily available tissues with varying light absorption properties. The three tissues chosen in this study have different light absorption and acoustic attenuation properties, but all achieved the expected results. Furthermore, we employed a laser diode system with an energy output of 1 $\mathrm{\mu}$J/pulse, producing significantly lower photoacoustic signals, in contrast to the typical high-power lasers used in prostate studies, where the energy output is typically measured in mJ/pulse, as shown in [34]. In future studies, we plan to investigate our low-energy laser diode for in vivo animal study. All three tissues yielded the anticipated results, showing the versatility of our methodology. This proof-of-concept study is essential for obtaining approvals for our future in vivo research using perfusion and oxygenation.

To evaluate the effect of poor signal on our search algorithms’ accuracy, we have performed an experimental analysis to find the correlation between the signal-to-noise ratios (SNR) and search errors. In PA imaging, averaging is the standard method to increase SNR. Therefore, we obtained different SNRs by controlling the number of average frames. Here, we defined SNR as [39]. As shown in Fig. 11, the search error is stable once the SNR exceeds 11 dB. In conclusion, the SNR will not much affect the search accuracy as long as it is larger than a minimum value where the PM starts appearing in the DAS images. In our experiments (online and offline), all PA images’ SNR exceeded 11 dB after 145 frames averaged, even though a much larger number of averages is achievable with our laser diode, enabling close-to-real-time performance at a repetition rate in the MHz range.

 

Fig. 11. The curve of SNR and search error. The raw data was collected on bovine muscle stained with ICG.

Download Full Size | PDF

5. Conclusion

In this study, we presented an automatic search method for PM to deal with the drawbacks of manual search procedures. We cast the PM search procedure as an optimization problem of a 1-D unimodal function and investigated three algorithms to address this problem. Simulation studies and experiments on various ex vivo tissues and phantom validated the validity of our proposed method. In simulation studies, all three algorithms perform well. In offline analysis, the average error of the TS and GSS algorithm on all materials was below 0.5$^\circ $ when acquiring 28 PA acquisitions. In online analysis, the average error of the TS and GSS algorithm achieved 0.29°${\pm} $0.31° and 0.48°${\pm} $0.32° after 28 acquisitions, while the WA achieved 0.53°${\pm} $0.30° using only 8.7 times acquisitions. In summary, the proposed technique significantly improved the system’s degree of automation and alleviated the labor required, while achieving good search accuracy. In future work, the proposed search method will be combined with real-time photoacoustic imaging to significantly improve the time efficiency of the PM-based RAPL.