1.

Introduction

Second-harmonic generation (SHG) imaging, based on the second-order nonlinear optical process, in which a noncentrosymmetric material (i.e., a material showing no center of inversion symmetry) converts a portion of incident light to scattered light at exactly twice the frequency, has developed into an important nonlinear imaging technique over the past several years. Among its benefits is its inherent ability to produce volumetric images of biological tissues comprising collagen fibers¹ without the need for labeling with an exogenous contrast agent. We have previously demonstrated quantitative SHG (Q-SHG) imaging as an effective and accurate modality to ascertain quantitative information from such collagen-based biological tissues.^2–6 For example, with respect to collagen fiber organization, the preferred orientation and orientation anisotropy have provided significant information.² Indeed, application of quantitative SHG has resulted in assessment of microstructural information of nonpregnant rat cervical tissue, healthy from injured horse tendons,⁷ age-related changes in porcine cortical bone,⁸ and dissimilarities in stromal collagen fiber organization in human breast biopsy tissues at various pathological stages.⁵ Additionally, Q-SHG imaging has also been reported as a suitable method for quantifying the change in dermal collagen fibers in skin burns using a rat skin burn model.⁹ In spite of these advancements, full applicability of Q-SHG imaging has been confined to static imaging conditions, thereby leaving its utility for dynamic biological processes largely unexplored. To address this, we recently reported on the experimental and computational requirements for carrying out Q-SHG imaging under dynamic conditions,¹⁰ i.e., simultaneously computing and displaying quantitative information with image acquisition. We found that for a $512\times 512\text{-}\text{pixel}$ area, the preferred orientation of collagen fibers in a tissue specimen captured by SHG imaging can be computed within $\sim 950\mathrm{ms}$ using a standard multicore CPU.

Recently, graphics processing units (GPUs) have emerged as alternate computation devices for faster processing compared to standard CPUs. GPUs comprise several thousand times more processing units or cores compared to conventional computer CPUs, permitting all cores to be used to carry out the same desired instructions in parallel. This facilitates its usage in various computational imaging applications where processing time is expensive. For example, GPU-based algorithms have been used to develop spectral (Fourier) domain optical coherence tomography^11–14 techniques with a significant reduction in computation time when compared to standard CPU-based algorithms. Similar results have also been observed from using the GPU for reconstruction of x-ray computed tomography images of high contrast and precision¹⁵ as well as performing deconvolution of three-dimensional confocal microscopic images.¹⁶ Additionally, GPUs have also been used to accelerate and optimize Monte Carlo simulations,^14,17 typically used to study the theory of light transport through various media. As such, the GPU would be extremely useful in obtaining quantitative information at the time scales of some of the faster biological processes, such as the propagation of an action potential in neurons occurring on the order of milliseconds, which has been successfully captured with SHG.¹⁸ The approach could also be useful for situations where analysis of multiple, quantitative metrics are incorporated with simultaneous image acquisition. In this work, we incorporate an NVIDIA GPU to our image analysis system, enabling parallel processing of our dynamic SHG image analysis algorithm. As proof-of-concept, we present several synthetic experiments in which we apply our GPU-based approach to quantitatively analyze consecutive frames from several videos of $512\times 512\text{-}\text{pixel}$ SHG images. In general, the videos are of varying arrangements of collagen fiber organization. We compare the computation time obtained using GPU versus CPU. The paper is organized as follows. Section 2 describes the experimental methods and image analysis technique used. Section 3 presents the results and discussion, while Sec. 4 provides the conclusion.

2.

Methods

2.1.

Image Analysis

We have previously provided a detailed description of our quantitative SHG image analysis carried out under dynamic conditions, which can be found elsewhere.¹⁰ Briefly, after acquiring a $512\times 512$ image, a Gaussian filter is applied and the image is subsequently divided into a $16\times 16$ grid, with each grid containing $32\times 32\text{pixels}$. For each grid, a preferred orientation is calculated based on the computed intensity gradient for each pixel within a grid. This information can then be used to estimate a global preferred orientation for the whole image. The accuracy of the calculated orientations is estimated by the circular variance,^19–22 a detailed description of which is provided in the Appendix. An intensity threshold is set to discriminate the background from the signal in a manner analogous to what we have previously reported.² This same threshold is used for calculations that are done both within a grid and for the global orientation estimate. Finally, an image is displayed with a gridded overlay, with arrows indicating preferred fiber orientation within each grid and the computed average orientation and circular variance being presented.

MATLAB® coupled with the compute unified device architecture (CUDA) parallel computing platform was used to develop the code. The parallel instructions were written in the C programming language, using the NVIDIA CUDA library version 5.0, and implemented in the GPU while MATLAB® was used as the host function to acquire the image, transfer image data to and from the GPU, and subsequently display the results. The hardware used for implementation consisted of an NVIDIA GTX 590 GPU, running on a Windows 7, Core i7-2600K Quad Core CPU running at 3.40 GHz clock speed, 3.8 GHz of maximum Turbo frequency, and 24 GB of DDR3-1066/1033 RAM. This same computer was used for the comparisons where GPU-based calculations were compared with CPU-based ones. A description of the GPU architecture and the CUDA programming model can be found in the CUDA Programming Guide 4.0.²³ To facilitate its adaptation in the GPU architecture, the image analysis procedure was divided into three segments known as CUDA kernels as shown in Fig. 1. In the first kernel, the acquired $512\times 512$ image is divided into a $43\times 43$ image grid, each of which is $12\times 12\text{pixels}$ in size. For this kernel, the GPU grid contains $43\times 43$ threadblocks, where each threadblock contains $16\times 16$ threads. Each threadblock is assigned to apply a Gaussian filter over one image grid. As the pixels in the boundary of a grid require contributions from neighboring pixels to apply the Gaussian filter, the threadblocks contain one extra thread in each boundary. In the second kernel, the filtered image is divided into a $16\times 16$ grid of $32\times 32\text{pixels}$ each, where the preferred orientation is calculated for each individual grid. In both these two kernels, individual GPU threadblocks are assigned to individual image grids, where individual pixels inside the grids are computed upon by individual threads in the threadblock; hence, parallel operation of data points is achieved on two separate levels. On the first level, all the threadblocks in the GPU multiprocessors carry out their computation in parallel, while on the second level, the individual threads inside a threadblock also operate in parallel. This essentially means that all pixels in a grid and all grids in the image are computed simultaneously, thus significantly reducing computation time. Finally, in the third kernel, the preferred orientations from the individual grids are used to calculate a global preferred orientation. As there are only 256 ($16\times 16$) preferred orientation values, a single GPU threadblock of 256 threads is sufficient. After performing the processing in the GPU, the image along with the quantitative information is returned to the CPU. MATLAB® instructions are used to display the preferred orientation in each block of the $16\times 16$ grid and a circular histogram showing the distribution of the preferred orientation values over the complete image. Finally, the global preferred orientation and the associated circular variance are also displayed.

Fig. 1

Flow chart of the steps performed in the graphics processing unit (GPU)-based quantitative image analysis modality.

3.

Results and Discussion

Figure 2 depicts the representative frames from the first video, as well as the results obtained by operating the two different computational modalities on it. Figure 2(a) shows the first six frames of the video, while Fig. 2(b) shows the frames that are captured and analyzed by the GPU-based and CPU-based codes in two consecutive rows, respectively. It is clear from Fig. 2(b) that the GPU-based code successfully captures and analyzes all six consecutive frames of the video. For the given video frame rate (10 fps), this is consistent with the expected computation time of $\sim 100\mathrm{ms}$ for each image. In comparison to this result, the CPU-based code only acquires the first frame and fails to capture any of the subsequent frames. From Video 1, it is also observed that the CPU is successful in capturing the first and the tenth frames of the video. This observation supports our previously reported fact that the computation time for a CPU-based code is $\sim 950\mathrm{ms}$ for a $512\times 512\text{-}\text{pixel}$ image,¹⁰ indicating that at 10 fps, it would fail to capture the next eight consecutive frames. The same analytical steps were carried out for the second sample video, the results of which are depicted in Fig. 2(c). The two consecutive rows in Fig. 2(c) show the frames that are captured by the GPU-based and CPU-based modalities, respectively. It is observed from Fig. 2(c) that for the new video frame rate of 33 fps, the GPU fails to capture two out of every three frames of the video. It is also observed from the second row in Fig. 2(c) that the CPU-based code could only analyze the first frame in this case, too, and fails to capture any subsequent frames shown in this figure. Video 2 reveals that the next frame successfully captured by the CPU is frame 32.

Fig. 2

Representative frames from a video of second-harmonic generation (SHG) images of breast biopsy tissues. (a) The first six consecutive frames in the video, showing fibers that are progressively rotated by 20 deg with respect to the horizontal. The frames captured and analyzed by the GPU-based and CPU-based image analysis are shown in (b) for 10 fps and in (c) for 29 fps. The scale bar corresponds to $10\mu \mathrm{m}$ for all images. These results are also compiled in the videos (Video 1 MPEG, 4.2 MB [URL: http://dx.doi.org/10.1117/1.JBO.19.9.096009.1] and Video 2 MPEG4, 7.3 MB [URL: http://dx.doi.org/10.1117/1.JBO.19.9.096009.2]).

Figure 3 shows the results obtained from analyzing consecutive frames of a video comprising SHG images of a variety of collagen-based tissues. Here, the goal is to evaluate the performance of the GPU-based code when SHG images vary (in fiber density, orientation, and organization) greatly between each consecutive image. Again, we observe in Fig. 3(b) that the GPU-based code captures and analyzes all consecutive frames from the video irrespective of fiber orientation or density, while in Fig. 3(c) and Video 3, we see that the CPU-based code could only perform its analysis on the first and the eleventh frames. Thus, the processing time for the GPU-based approach is not influenced by fiber density and spatial organization.

Fig. 3

Representative frames from a video of SHG images of several collagen-based tissues (a). Frames 1, 3, and 4 are SHG images of human breast biopsy tissues, while frames 2 and 5 are of rat cervix and porcine tendon tissues, respectively. (b) and (c) show the frames captured and analyzed by the GPU-based and CPU-based codes, respectively. These results are also compiled in a video (Video 3 MPEG4, 5.3 MB) [URL: http://dx.doi.org/10.1117/1.JBO.19.9.096009.3]. The scale bar corresponds to $10\mu \mathrm{m}$ for all images.

Figure 4 compares the performance in processing time between the GPU- and CPU-based approaches for each segment of the analysis. To get an accurate estimate of the processing times, the two modalities are used to analyze 20 SHG images of size $512\times 512\text{pixels}$ showing varying degrees of fiber organization and density. The processing times for each step for all the images are obtained and their average values are used in this comparison. It is observed from Fig. 4 that the calculation of the preferred orientation of individual image grids (kernel 2) takes the longest amount of time for both approaches. As such, to obtain a significant reduction in overall computation time, it would be important to reduce the time required by kernel 2. Our GPU-based implementation achieves a time improvement of $\sim 20\times$ for kernel 2, reducing it from $\sim 520\mathrm{ms}$ used by the CPU to $\sim 18\mathrm{ms}$. Application of a Gaussian filter and calculation of the global preferred orientation were performed $\sim 5\times$ and $\sim 35\times$ faster, respectively. The time required to display the results (not shown) is the same for both the modalities and it was observed to be $\sim 50\mathrm{ms}$. Overall, the GPU-based code performs the analysis at an average of $\sim 10\times$ faster than the CPU-based code.

Fig. 4

Image processing time required for the GPU-based analysis compared with the CPU-based analysis for three different analysis steps.

It is worth noting that although $512\times 512\text{-}\text{pixel}$ images were used for the proof-of-concept, the GPU-based code can also be used to analyze images of higher pixel density without any modification. The image grid size and the number of threadblocks in the GPU would scale according to the image size. Individual image grids and threadblocks would also contain the same number of pixels and threads, respectively. However, this would not be possible in the current version of the GPU due to limitations in the number of threadblocks and grids that can operate in parallel. In certain cases, the individual image grid sizes may need to be increased to facilitate clear visualization. This would require each thread in a threadblock to process more than one pixel. For example, if the image grid contains $64\times 64\text{pixels}$, a threadblock of $32\times 32$ threads would assign one thread to analyze data obtained from four pixels. Note that this would require an increased amount of memory in the GPU, which is not supported in the current version of the GPU that was used in this work. Apart from this, if any further quantitative analysis is desired from the SHG images, additional CUDA kernels can be constructed and conveniently added to the existing code.

4.

Conclusion

In this paper, we demonstrated GPU-based quantitative analysis of SHG images. We showed that the preferred orientation of collagen fibers can be determined in $\sim 100\mathrm{ms}$, either at the level of individual elements in a $16\times 16$ grid or globally for a $512\times 512\text{-}\text{pixel}$ image. As proof-of-concept, we have applied this modality to analyze consecutive frames of two videos of 10 and 33 fps, respectively. In the first case, the GPU-based system successfully captured and analyzed all the frames of the video, while in the second case, it succeeded in capturing one in every three frames. In contrast, the same analysis using a standard CPU failed to capture all but the first frame from both the videos for the representative frames shown. Both approaches were again compared for a video of SHG images from more complex collagen-based structures, with the GPU-based method clearly outperforming the standard method. This improvement in processing time makes our approach attractive compared to other nonlinear imaging modalities that are modified for quantitative analysis.