Tattoo tomography: Freehand 3D photoacoustic image reconstruction with an optical pattern

Purpose Photoacoustic tomography (PAT) is a novel imaging technique that can spatially resolve both morphological and functional tissue properties, such as vessel topology and tissue oxygenation. While this capacity makes PAT a promising modality for the diagnosis, treatment, and follow-up of various diseases, a current drawback is the limited field of view provided by the conventionally applied 2D probes. Methods In this paper, we present a novel approach to 3D reconstruction of PAT data (Tattoo tomography) that does not require an external tracking system and can smoothly be integrated into clinical workflows. It is based on an optical pattern placed on the region of interest prior to image acquisition. This pattern is designed in a way that a single tomographic image of it enables the recovery of the probe pose relative to the coordinate system of the pattern, which serves as a global coordinate system for image compounding. Results To investigate the feasibility of Tattoo tomography, we assessed the quality of 3D image reconstruction with experimental phantom data and in vivo forearm data. The results obtained with our prototype indicate that the Tattoo method enables the accurate and precise 3D reconstruction of PAT data and may be better suited for this task than the baseline method using optical tracking. Conclusions In contrast to previous approaches to 3D ultrasound (US) or PAT reconstruction, the Tattoo approach neither requires complex external hardware nor training data acquired for a specific application. It could thus become a valuable tool for clinical freehand PAT.


Introduction
Photoacoustic tomography (PAT) is an emerging imaging modality that has undergone rapid development since the late 90s [1]. PAT is based on the photoacoustic (PA) effect [2] and, using a light in-sound out principle, provides high contrast in optically absorbing tissue regions. The tissue is illuminated with nanosecond-long pulsed near-infrared laser light. When the light is absorbed by specific molecules, thermoelastic expansion causes a local pressure rise that generates an acoustic shock wave which can be measured by a standard ultrasound (US) transducer. Furthermore, acquiring multispectral PA data allows not only to spatially resolve highly concentrated absorbing tissue regions, but also to define the absorbing molecules and estimate their concentration in a process referred to as spectral unmixing [3]. PAT thus enables the measurement of both morphological and functional tis-sue information in depths up to several centimeters. A key application of PAT is the estimation of blood oxygenation, which is relevant for a variety of diseases [4].
Clinical PAT devices typically use handheld 2D ultrasound transducers for data acquisition. The bottleneck of such systems is that they produce only cross-sectional 2D slices of the imaged structures, although the full 3D context is crucial for various applications [5,6]. Current hardware solutions toward 3D imaging include whole-body scanners-which only allow small animal imaging-or 3D probes, which lack penetration depth, field of view, and spatial resolution and are associated with additional costs [7]. Another approach toward 3D US/PAT is based on the constrained movement of a 2D probe, either introduced by additional mechanical hardware [8] or lightweight robots [9], which also offers additional intrinsic tracking data. Software-based solutions are an alternative for compounding 3D volumes based on pose estimation of the individual 2D image slices. However, these may also rely on complex hardware, such as an external tracking system [10]. Mobile cameras mounted to the probe for motion tracking have also been investigated [11] but come at the cost of a complex hardware setup as well. To overcome these issues, sensorless approaches based on speckle decorrelation [12,13] have been proposed for US, but they are not well applicable to PA images, which do not show speckle effects due to the physical difference in image formation between PAT and US [14]. Recent work enhanced the speckle-based approach with deep learning-based components [15]. However, one of the core limitations of machine learning-based approaches is that they are commonly trained for a specific application using a tailored training data set [15,16]. In general, this leads to a lack of robustness and generalizability and, therefore, poor performance when transferred to new settings [17,18]. While this issue may be addressed in future work, the purpose of the present work was to address current problems with an entirely novel approach that specifically leverages the characteristics of PA imaging. It should be pointed out that we are not aware of any prior work in sensorless freehand 2D to 3D PAT reconstruction, while the reader may refer to previous work in freehand 3D US reconstruction [13,15,19] for a broader review on the topic of 3D image compounding.
To overcome the limitations of prior work, we present an entirely novel approach to 3D reconstruction that leverages the unique capability of PA to tomographically reconstruct optical properties. Specifically, the concept is based on an optical pattern placed on the region of interest which allows the recovery of the 3D probe pose relative to the pattern coordinate system based on a single tomographic image. By transforming a sequence of freehand PAT images into the common pattern coordinate system, a consistent 3D reconstruction of PAT imaging data can be obtained without the need for an external tracking system.
The following sections present the prototype implementation of the novel approach (patent pending) along with a feasibility study assessing the accuracy and precision of the reconstruction.

Material and methods
In this section, (1) the general concept of Tattoo tomography is outlined, followed by a description of the prototype implementation, comprising information on (2) the used PAT imaging device, (3) the first optical pattern designed for the approach, (4) the corresponding approach to pose estimation, and (5) the method for compounding. Finally, (6) the validation concept is described.

Tattoo tomography concept
The process of reconstructing a 3D image with the Tattoo tomography concept is illustrated in Fig. 1. 1. Pattern attachment: Prior to image acquisition, an optical pattern is placed above the anatomical region of interest. This pattern should fulfill the following key requirements: (1) It should be designed in a way that a tomographic image of it enables the estimation of the probe pose relative to the coordinate system of the pattern and that (2) it is easy to mount and remove in a clinical setting. (3) The dye that makes up the pattern (here seen in green) should absorb in the frequencies that match those of the imaging system (here: near-infrared) and that are complementary to those used for the actual imaging. 2. Image acquisition: A sequence of images is acquired, each showing a part of the pattern as well as the target region. 3. Pose estimation: As the pose information is encoded in the optical pattern by design, the pose of each image slice within the pattern coordinate system can be determined via the PA image features corresponding to the pattern. 4. Image compounding: The 3D volume is reconstructed by interpolating between the transformed slices in the pattern coordinate system.
The following sections present the first prototypical implementation of this concept.

PA imaging device
All PA images used as part of this work were acquired with a multispectral optoacoustic tomography (MSOT) Acuity Echo research system (iThera Medical GmbH, Munich, Germany). It is a hybrid device, which can acquire co-registered and synchronized PA and US images. The built-in laser is (2) Next, a sequence of images is acquired, each showing a part of the pattern as well as the target region. (3) The image features corresponding to the pattern (here: three points) are extracted from each slice and used to recover the pose of the probe relative to the pattern coordinate system. (4) Once all acquired slices have been transformed (T) into a common coordinate system (C), the 3D volume is compounded an Nd:YAG laser with a tuning range from 660-1300 nm, a peak pulse energy of 30 mJ, a repetition rate of 25 Hz, and a pulse duration of 4-10 ns. The concave 2D US transducer has a 4 MHz center frequency, 256 elements, and a radius of 40 mm.

Optical pattern
For this first prototypical implementation, a trident design was chosen, as depicted in Fig. 2. This trident consists of an isosceles triangle (for simplicity referred to as tilted lines) and the corresponding angle bisector between the tilted lines of the triangle (referred to as central line).
Mathematically speaking, this initial prototype pattern design currently only allows accounting for three degrees of freedom (DoF), representing the intersection line of the image coordinate system with the tattoo coordinate system. The remaining three DoFs are recovered by imposing constraints on the image acquisition process: Firstly, the optical pattern has to be approximately placed such that it resembles a flat surface, and secondly, the PA probe is required to be positioned orthogonally to the pattern plane (orthogonality constraint).
The pattern was printed on a transparent foil commonly used in overhead projection (PrintLine Overheadfolie, Vleveka, Deurne, Netherlands) using cyan ink of a print station (TASKalfa 5052ci, Kyocera, Esslingen, Germany). These foils are well suited for skin contact, easily sterilizable, and printable in a straightforward manner. During acquisition, the pattern provides a strong contrast at a wavelength of 750nm, while the functional imaging was performed at a wavelength of 850nm, where oxygenated hemoglobin is a dominant absorber in tissue. The foil is placed on tissue with a thin layer of coupling gel on top and underneath the foil and fixated using tape to avoid movement of the foil relative to the skin. The setup allows to easily place and remove the pattern on/from the tissue of interest.

Pose estimation
The pose estimation step refers to estimating the pose of an individual image slice relative to the pattern coordinate system (C tattoo ). To this end, the information corresponding to the optical pattern is extracted from the PAT image that reflects the wavelength in which the pattern is absorbing (here: 750 nm). When using the optical pattern proposed in the previous section, this task is equivalent to extracting the three intersection points of the image with the line-based pattern, as depicted in Fig. 1. This is achieved with the SciPy [20] function find_peaks. Based on the three points, the line corresponding to the probe pose on the optical pattern (blue line in Fig. 2) can be unambiguously determined. It is characterized by the distance a 0 of the absorption peak of the central pattern line c to the pattern origin and the slice tilt angle α as defined in Fig. 2. These are computed as follows: where d r and d l are the distances between the right and left point to the central point, respectively, and γ is the pattern opening angle. When assuming an orthogonal pose of the probe relative to the pattern plane, the line given by α and a 0 uniquely determines the pose of the probe in 3D. This results in the following transformation matrix T img− →tattoo which transforms the PA image slice into the tattoo coordinate system: where x c and y c are the x and y coordinates of the central tattoo point c in the PA image slice and α and a 0 the previously derived tattoo coordinates.

3D volume compounding
The 3D compounding algorithm is inspired by the "Volume reconstruction algorithm" of the public software library for ultrasound imaging research (PLUS) [21]. A 3D zero-valued voxel grid is defined inside the bounding box of all transformed slices. The voxel spacing is set corresponding to the pixel spacing of the image in x-y-direction and predefined according to the probe movement in z-direction. Each image pixel in each slice is assigned to the closest voxel location. The individual voxel values are then calculated as the average pixel value of all assigned pixels.

Experimental setup
The purpose of our experiments was to assess the general feasibility of the Tattoo tomography concept. We provide (1) a qualitative visual demonstration, (2) a quantitative validation of the accuracy in a phantom experiment, and (3) a quantitative validation of the precision in an in vivo experiment. The pattern was constructed with an opening distance d max of 20 mm and a vertical extension (a 0,max ) of 50 mm. The scans of all experiments were performed using the two wavelengths 750 nm and 850 nm, where 750 nm corresponds to the wavelength in which the optical pattern is absorbing strongly (and thus well visible) and 850 nm corresponds to the wavelength used for the construction of the functional tissue information. The PA data were reconstructed using a delay-and-sum [22] algorithm.
Feasibility demonstration For a demonstration of the principal feasibility of the approach, a logo of the international conference on Information Processing in Computer Assisted Interventions (IPCAI) was printed and placed inside a gel pad (AQUAFLEX Ultrasound Gel Pad, Parker Laboratories, inc., Fairfield, USA) that was cut in half in height (see Fig. 4 a,b). The optical pattern was fixed on top of the gel pad with pins. For scanning, the probe was moved along the pattern in two principal ways: (1) slow and without direction changes, representing a careful acquisition and (2) fast with many direction changes, representing an extremely careless acquisition (see Fig. 5, bottom row). It has to be noted that, performing the freehand scans (for the feasibility and in vivo study), the probe was only approximately orthogonal. This means that the orthogonality of the probe was dependent on the skill of the scanning person and might be distorted to a small degree reflecting a realistic setup.
Phantom validation The phantom used for the phantom validation was an N-wire phantom (see Fig. 3) typically used for hand-eye calibration of tracked devices [10,21]. The phantom was built using three layers of N-shaped black nylon wires (0.4 mm in diameter) showing a strong broadband absorption in the near-infrared. It was placed within a milk bath containing 1.5 l of tap water and 20 g of coffee cream (10 % fat) to increase the scattering coefficient at room temperature. The phantom wiring was transferred to a 3D model based on the known location of the wires. A point representation of this model was then used as a reference for the validation, as detailed below (see Fig. 6).
For image acquisition, the optical pattern was fixed on top of the phantom using Velcro ® strips. This setup was placed upon a movable optical breadboard, and the probe was clamped above the phantom with a fixed probe angle α (rotation around y-axis). Sensorless approaches to 3D US reconstruction typically rely on the presence of speckle effects. As PA images do not feature these effects [14], US reconstruction methods are not directly applicable to PAT data. For this reason, we compare the Tattoo approach to a reference method based on optical tracking, which is the only freehand PAT method in the field of 2D to 3D reconstruction that we are aware of [10] and which is a widely used baseline method in 3D US reconstruction [15] [19]. As proposed in [10], the probe was tracked with an optical tracking system (NDI Polaris Spectra camera, NDI Medical, Waterloo, Canada) and calibrated accordingly, where a marker was attached to the probe and a stationary marker was placed on a table next to the setup defining the world coordinate system. The temporal calibration was performed via registration of image and tracking data of a sinusoidal scan of the ground of a box filled with water, allowing to temporally match the movement of the probe (tracking data) and the distance of images represent an intersection of the image stack at a fixed y-position before (left) and of the 3D volume after (center) Tattoo reconstruction Fig. 6 The 3D model of the N-wire phantom (white) is displayed together with the optical pattern compounded volumes of a single scan (blue, left) and all three scans (blue, red, green, right) after applying an ICP algorithm to register the Tattoo point cloud on the N-wire model (cf. Table 1) the probe to the ground (image data) [23]. The spatial calibration was accomplished using an N-wire phantom [21], deriving the transformation between markers attached to the probe and the image plane by registration of known phantom wire position [24]. During the scanning, the breadboard was moved under the probe to allow for a smooth acquisition of the area of interest. The PAT scans were then compounded using both the Tattoo method (cf. Sect. 2.5) and the optical tracking baseline method described in [10]. Subsequently, a 3D segmentation of the PA wire volume was performed by applying a background suppressing threshold using the Medical Imaging Interaction Toolkit [25].
Based on the known N-wire model and the 3D reconstructions, the accuracy of the compounded images was assessed. As the relation of the phantom coordinate system and the tattoo coordinate system can only be approximated, we concentrated the validation on the consistency of the reconstructed model. Specifically, we applied the iterative closest point (ICP) algorithm as implemented in the MITK [26] framework (see Fig. 6) to register the reconstructed point cloud with a (simulated) point cloud representing the N-wire model. The fiducial registration error (FRE) was then used for quantifying the accuracy of the 3D reconstruction (see Table 1). As the stopping criterion for the ICP algorithm, the change of the FRE between two ICP iterations had to be smaller than 0.001 mm.
In vivo validation In the in vivo validation step, we concentrated on validating the precision (reproducibility) of our method. To this end, we obtained a total of 30 tattoo scans from three healthy volunteers. For each volunteer, the optical pattern was placed on the forearm, and 10 measurements were taken while tracking the PA probe with an optical tracking system. All measurements were transformed to the tattoo coordinate system (for our method) and the coordinate system of the optical tracking system (for the baseline method), and two corresponding 3D reconstructions were performed. The largest image region, corresponding to a vessel and comprising measurements of at least seven volumes, served as the target region for the assessment. To approximate precision, we performed a slice-based analysis. The reference position of the vessel center was computed as the mean of the (up to) ten vessel locations extracted from the individual scans. The latter were again determined with the find_peaks function (see Sect. 2.4). The mean Euclidean distance between this reference and the individual measurement served as the metric for the precision of the compounding and referred to as mean vessel distance (MVD).  Figures 4 and 5 show the setup and the results for the feasibility experiment. It can be seen that the Tattoo reconstruction is feasible even when images are acquired in an extremely careless and unsystematic manner.

Phantom validation
The results of the Tattoo-compounded volumes registered with a 3D model of the N-wire phantom can be seen in Fig. 6. The FREs for the registration of the compounded volumes (Tattoo and optical tracking) with the N-wire model are shown in Table 1 for three different, but fixed, rotation angles α (0 • ,4 • ,8.5 • ). It is 0.63 mm on average for the Tattoo compounding and 0.87 mm on average for the optical tracking-based compounding.
In vivo validation For one volunteer, it was not possible to clearly identify a continuous vessel within at least 7 out of 10 3D volumes to satisfy our metric condition (see section 2.6 )-neither in the baseline volumes nor in the Tattoo tomography volumes. Therefore, our results are restricted to the 20 scans of the remaining two volunteers, referred to as v1 and v2. As summarized in Table 2, the Tattoo-based registration is slightly more reproducible compared to the optical tracking-based method.

Discussion
In this paper, we presented Tattoo tomography, a novel approach to 3D PA image reconstruction from an acquired sequence of 2D PA image slices. The key advantage of the concept is its simplicity; applying a pattern without the need for involved calibration, training of an algorithm, or any additional complex hardware allows smooth integration of the proposed approach into any clinical workflow. Our prototype implementation of the proposed concept further illustrates that such a pattern-based concept may also feature higher accuracy and precision compared to widely used techniques that rely on external tracking devices. An important advantage of the Tattoo approach is that the clinically relevant tissue information and the information needed to decode the probe pose can be encoded in different wavelengths. We thus do not expect reduced image quality resulting from the attachment of the pattern. It deserves mentioning in this context that we chose the materials (ink and foil) of the pattern such that the influence on acoustic and optical propagation within the functional PA wavelength range is kept minimal and we did not qualitatively observe any reduced image quality.
It is worth mentioning that we had initially planned to use the optical tracking-based approach as a reference method. First experiments, however, suggested that the proposed approach might be even more accurate and precise than this method, which is why we decided on a comparative assessment. The FRE obtained for both reconstruction methods may appear relatively high (0.6-0.8 mm). We attribute this to the fact that the model representation of the phantom was only an approximation of reality. This hypothesis is supported by the fact that both reconstruction methods feature a systematic mismatch of the outer wires, as shown in Fig. 6.
In our in vivo experiment, we examined the precision of Tattoo tomography in comparison with the optical tracking baseline. The fact that the Tattoo approach was substantially more precise than the baseline method in terms of the MVD may most likely be attributed to the different underlying world coordinate systems. The tattoo coordinate system is a relative coordinate system that moves with the target structure and is, therefore, influenced by patient motion to a lesser degree. Due to the limited number of volunteers and the lack of a reliable reference for the reconstruction, further studies are required to support these initial results especially regarding the robustness, repeatability, and accuracy in vivo. Moreover, this then enables us to evaluate the approach within a specific clinical setting.
While the general Tattoo concept proposed is potentially very powerful, our first prototype implementation comes with several limitations. The main current drawback is the fact that we use a 2D pattern. Future work could be directed to developing a 3D pattern, allowing for the recovery of 6 DoF from a single measurement, thus eliminating the need for holding the probe orthogonally to the pattern. Furthermore, the trident design was chosen for its relatively simple optical pattern. While-in theory-it enables us to uniquely recover the pose of the device from a single measurement and thus fulfills the core requirement, more work needs to be invested in patterns that are optimized for the robustness of pose estimation in the presence of uncertainties (e.g., uncertainties with respect to point localization). In this context, dyes that absorb at different wavelengths of light may play an important role in the future.
Assessment of the flat surface constraint In this first feasibility study, we have assumed the optical pattern to be a flat surface. It should be noted, however, that violation of this constraint would potentially lead to errors in the pose estimation. The primary effect of a non-flat surface is the change of distances of the pattern absorption points, as illustrated in Fig. 7 (left). From a mathematical point of view, the flat surface assumption allows us to compute the Euclidean distance between points, while a curved surface would require more advanced methods, such as the computation of the geodesic distance along the skin surface.
To analyze the potential effect of violations of the flat surface assumption, we selected a PA image of our test set that featured a particularly high curvature of the tattoo. We then approximated the geodesic distance between the tattoo points in one image slice with a manually placed Bézier curve using two control points, as shown in Fig. 7. The resulting increase in distance compared to the Euclidean baseline was 2% (17.7mm vs. 17.4mm). This discrepancy corresponds to a shift of a 0 = 0.4 mm in the pattern coordinate system and an angular error of α = 1.7 • . This can be regarded as tolerable given the overall resolution of the imaging system and the fact that we picked one of the most extreme cases.
We see primary ways to address the issues arising from the flat surface constraint in the future: 1. Image processing-based solution: Automatically segmenting the skin from the image would enable us to approximate the geodesic distance along the skin surface between any two given tattoo points by determining the length of the curved line representing the skin in a 2D tomographic image. As isometric transformations such as the bending of a flexible material (here: the tattoo) preserve geodesic distances [27], a curved line represented by three tattoo points can be related to a corresponding straight line in the flat tattoo coordinate system via the Euclidean distance (which is equal to the geodesic distance on surfaces with zero curvature). Leveraging this principle, the curved tattoo could be reconstructed in 3D and serve as a reference coordinate system just like the flat variant. While this approach could compensate for small deviations from the flat surface constraint, a highly curved target anatomy would still pose a challenge when used in conjunction with a 2D pattern. In this case, the 3D reconstruction of the tattoo itself would be possible in theory, yet, compounding errors of the 3D anatomy could still occur as the transfer of the orthogonality requirement to a curved pattern is not straightforward. In this case, a hardware-based solution may help. 2. Hardware solution: Alternatively, a semirigid pattern could be applied, comprising a flat tattoo at the top and Finally, we are planning on enhancing the optical pattern, such that it can not only be used for 3D reconstruction but also for multimodal registration with conventional imaging data, such as computed tomography (CT) or magnetic resonance imaging (MRI) data.
In conclusion, we have presented the first 3D PAT reconstruction concept that elegantly leverages the unique capabilities of PAT. The proposed concept is simple to apply, does not require additional devices such as a tracking system, and can easily be integrated into clinical workflows. According to initial results, it may even be more accurate and precise than competing methods that rely on complex hardware setups. We, therefore, assume that the concept has a high potential for future clinical translation. Availability of data/code and material Within the current paper version, interested readers are encouraged to get in contact with the corresponding author in case further insight into methods, image data or evaluation code is desired.

Conflicts of interest
The authors declare that they have no conflict of interest.

Ethics approval
The healthy human volunteer experiments were approved by the ethics committee of the medical faculty of Heidelberg University under reference number S-451/2020, and the study is registered with the German Clinical Trials Register under reference number DRKS00023205.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.