Extendable, large-field multi-modal optical imaging system for measuring tissue hemodynamics.

Simultaneous imaging of multiple hemodynamic parameters helps to evaluate the physiological and pathological status of biological tissue. To achieve multimodal hemodynamics imaging with a large field of view, an infinite conjugate relay lens system compatible with the standard C-Mount camera lens is designed to adapt one camera lens with multiple CCD/CMOS cameras for simultaneously multi-wavelength imaging. Using this relay lens system, dual wavelength reflectance imaging and laser speckle contrast imaging were combined to simultaneously detect the changes in blood flow, oxygenation, and hemoglobin concentrations. To improve the accuracy of hemoglobin concentration measurement with an LED illumination source, an integral algorithm is proposed that accounts for the dependence of differential pathlength factors (DPF) on hemoglobin concentrations and the integral effect of both the emission spectrum of the light source and the spectrum response of the detector. The imaging system is validated by both phantom and in vivo experiments, including the arterial occlusion, and the detection of blood volume pulse (BVP) and blood flow pulse (BFP) signal in human subjects. The system helps in the exploration of macroscopic tissue hemodynamics.


Introduction
Functional optical imaging such as Laser speckle contrast imaging (LSCI) [1][2][3][4], optical intrinsic signal (OIS) imaging [5][6][7][8][9], and fluorescence imaging [10], which have the advantages of real time, wide-field and non-contact, have been widely attempted in clinic applications. Switching the illuminating sources in time sequence is widely used to realize multimodal imaging [11][12][13][14][15][16]. This method sacrifices temporal resolution so that the rapidly changing signals is difficult to be resolved. Another approach utilizes the multiple charge coupled device (CCD) or complementary metal oxide semiconductor (CMOS) imaging sensors [17][18][19][20][21][22]. Multi-modal imaging systems using multiple cameras are usually constructed based on microscopes, which is relatively easy to realize because of the extendable optical structure of the infinite conjugate microscopes. Since the parallel beam paths are adopted between the objective and the tube lens, beam splitting elements or filters can be easily inserted into, and the distance between the objective and tube lens can be adjusted as needed. For examples, Liu et al reported a head-mounted multi-modal imaging system to observe cortical hemodynamic in freely moving animals by employing an infinite conjugate microscope as the relay lens unit [17]. Cha et al reported a simultaneous polarimetry, laser speckle contrast and fluorescence imaging system for intraoperative visualization of peripheral nerves and micro-vasculatures based on a surgical microscope [18].
However, in clinical diagnosis and therapies, besides the hemodynamics changes in microvascular, the spatiotemporal changes of hemodynamics in large field of view is also important in many applications such as evaluation of burns [23], nerve block [24], and feedback of laser therapy [25]. Most commercialized camera lenses available with large field of view have standard interfaces like C-Mount, F-Mount, etc. The back focal length of these lenses are usually 12-40mm, which is usually not long enough to accommodate beam splitting elements or filters. Therefore, it is necessary to introduce an extra relay lens unit for multi-modality imaging systems using the strategy of multiple sensors. Themelis et al extended the back focal length of camera lenses by utilizing a finite conjugate relay lens system [22]. However, the relay distance of the finite conjugate relay lens system is limited by its focal length. Furthermore, the large aperture angle of beams in the finite conjugate relay lens system may affect the transmittance properties of some optical filters.
Simultaneously monitoring tissue hemodynamic parameters such as blood flow, blood oxygenation and blood volume can help to better differentiate the physiological and pathological status of tissues, which has important clinical significance [26]. Based on Lambert-Beer's law, accurate measurements of blood oxygen and hemoglobin concentrations depend on the estimation of molar extinction coefficients and differential pathlength factors (DPF) [27]. Previous study usually utilized monochromatic illumination sources, such as multiple lasers [15,28], or a combination of a broad-spectrum illumination sources and multiple filters [11,17]. Light-emitting diodes (LEDs) have the characteristics of low price, high brightness and biosafety, which are suitable for functional imaging. Due to the wide illumination spectrum of LEDs, the differences in molar extinction coefficients and DPF cannot be ignored within this range. Hashimoto and Kurata utilized modified extinction coefficients called average extinction coefficients (AECs) to correct the influence from the wide spectrum of illumination sources and image sensors [26,[29][30][31]. Wang et al estimated the DPF based on average hemoglobin concentrations in tissue, but only a small variation in hemoglobin concentrations was considered [28]. The hemoglobin concentrations of different human tissues vary greatly, approximately range from 0-117uM [26,32]. Therefore, the influence of hemoglobin concentrations should also be considered when estimating DPF.
In this study, to realize large-field of view multi-modality imaging of tissue hemodynamics, we designed an infinite conjugate relay lens system compatible with the standard C-Mount camera lenses, which extended the back focal length of the camera lenses so that the camera lens could operate with multiple cameras simultaneously and the relay distance could be flexibly adjusted. Besides, we introduced an electrically tunable focusing lens for rapid focusing, which was convenient for the practical applications. With this relay lens system, an RGB CCD camera coupled with a dual wavelength LED were used for hemoglobin concentrations and blood oxygen imaging, and combined with the laser speckle contrast imaging recorded by another monochrome camera. The system allowed simultaneous measurement of relative blood flow (rBF), oxyhemoglobin (HbO) concentration, deoxyhemoglobin (HbR) concentration, total hemoglobin (HbT) concentration and relative metabolic rate of oxygen (rMRO 2 ). We proposed an integral algorithm that considers the dependence of DPF on the hemoglobin concentrations and the integral effect of both the emission spectrum of light source and the spectrum response of detector so as to improve the accuracy of the estimation of hemoglobin concentrations.

Optical imaging system
In order to fully relay the image from the image plane of the camera lens to the sensor, a 1:1 relay lens system compatible with C-Mount camera lenses was designed with Zemax software. As shown in Figure1(a), the relay lens system consists of two set of relay lens units. Each relay lens unit consists of a Steinheil triplet achromatic lens (f = 50 mm Edmund Optics, United States) and a meniscus lens (f = -150 mm Thorlabs, United States). The meniscus lens is used to reduce the field curvature of the system. Two sets of the relay lens units were arranged symmetrically to form a 1:1 relay lens system with parallel beam paths between them. To enable the electrical focusing, an electrically tunable focusing lens (EL-16-40-TC Optotune, Switzerland) was added in front of the relay lens system.
The multi-modal imaging system using the above relay lens system is shown in Fig. 1(b). To achieve both a large field of view and a compact structure, a camera lens with a focal length of 8 mm (M0824-MP2, Computar, Japan) was employed. The sample was firstly imaged through the camera lens, and then transferred through the relay lens system to the imaging sensor. The illumination source composed of a near infrared laser diode (785nm 90mW L785P090 Thorlabs, United States) and a high-power color LED (XML-COLOR Cree, United States) which integrated red, green, blue, and white LED chips in one package. The LED beam and the laser beam were combined through a long pass dichroic mirror (650nm DMLP650 Thorlabs, United States) and then passed through an engineered diffuser (ED1-C50 Thorlabs, United States) to produce a light spot large enough to cover the field of view. Between the two sets of the relay lens units, a long pass dichroic mirror (650nm DMLP650 Thorlabs, United States) was used to separates the visible light from the near infrared light. A laser filter (Edmund Optics, United States) was placed in front of the monochrome camera CCD1 (2/3 inch, QE, PCO Computer Optics, Germany) which received near infrared laser speckle pattern to calculate blood flow. The color camera CCD2 (1/2 inch, avA1000-100gc Basler, Germany) was employed to receive OIS signals to calculate hemoglobin concentrations and blood oxygenation.
The working distance of the image system, which is the distance between the object plane and the camera lens, can be adjusted by changing the operating current of the electrically tunable focusing lens at the range of 90-280 mm, obtained by Zemax simulation and shown in Fig. 1(c). The angle of the field of view in the imaging system was mainly depended on the camera lens, which was approximately 52 degrees.
The sample was illuminated by the red LED (630 nm) and green LED (530 nm) simultaneously. The CCD sensor in the color camera was equipped with a Bayer filter, which has red, green, blue channels, as shown in Fig. 1(d). The dual wavelength OIS signals were separated by the Bayer filters. To reduce the crosstalk between the OIS signals of two wavelengths, the images obtained by the blue channel and red channel of CCD2 were extracted to calculate the hemoglobin concentrations and blood oxygenation, because the spectrum responses of these two channels were less aliased.
The crosstalk was tested by using a white diffusive reflectance standard (Chuangpu, China) as well as the same red LED and green LED used in the following experiment. The ratio of the CCD counts of red channel to that of blue channel was 17.55 when the reflectance standard was illuminated by the red LED. The ratio of the CCD counts of blue channel to that of red channel was 9.02 when the reflectance standard was illuminated by the green LED.

Data acquisition and processing
A control software based on Labview (National Instrument, the United States) was programmed to trigger the cameras and operate the electronically tunable focusing lens. A custom-written Matlab (MathWorks, the United States) code was employed to register the images and calculate the hemodynamic parameters.
The data processing procedure is shown in Fig. 2. Blood flow was obtained by laser speckle contrast analysis using the speckle images captured with the monochrome camera CCD1. The blood oxygenation and hemoglobin concentrations were calculated with OIS images captured by the color camera CCD2. Then rMRO 2 was calculated using the results of rBF and hemoglobin concentrations after the images obtained by two cameras were registered.

Hemoglobin concentrations calculation
For a semi-infinite medium, when photons are assumed to come from an instantaneous point source at a depth z 0 below the surface, the DPF D a can be expressed as follows [26,33]: where u a is the absorption coefficient and u s is the reduced scattering coefficient of the medium. Most biological tissues are characterized by strong optical scattering [34]. When u s >> u a , it can be considered that: The relationship between u a and hemoglobin concentrations is shown as below: If the influence of oxygen content on D a was ignored, the relationship between D a and hemoglobin concentration can then be expressed as: c HbT = c HbO + c HbR (6) where α is a scale coefficient.
The changes in the reflected light intensity are related with the changes in hemoglobin concentrations by the modified Lambert-Beer's law [27]: where R 0 is the reflected light intensity of the baseline, R t is the reflected light intensity at time t.
The ε HbO and ε HbR are the molar extinction coefficients of HbO and HbR at a specific wavelength, respectively. Here, the DPF D a was determined by optical properties of tissues [32] using Monte Carlo simulation [35]. The CCD counts measured by a channel of the color camera can be expressed as [36]: where r is the spectrum response of the channel, R is the reflected intensity, L is the spectrum distribution of the illuminating light source and m is the analog to the digital conversion factor of the camera.
Considering the spectrum distribution of the illuminating source and the cut-off wavelength of the dichroic mirror, the integration was carried out in the range of 400-650nm. Therefore, using the first order approximation of Taylor expansion of [Eq. (7)], the calculation formula of hemoglobin concentration based on the color camera detection and dual wavelength LEDs illumination can be approximately expressed as follows: ∫ 650 400 r B L g ε HbR αdλ ∫ 650 400 r R L r ε HbO αdλ ∫ 650 400 r R L r ε HbR αdλ where, V B,0 , V B,t represent CCD counts of the blue channel at the baseline and time t, respectively. V R,0 , V R,t represent intensity of the R channel at the baseline and time t, respectively. m is the analog to the digital conversion factor of the camera that can be estimated by baseline conditions. r B , r R are the spectrum responses of the blue channel and the red channel respectively. L g , L r are the spectrum distribution of the green LED and the red LED. The constant α was obtained by Monte Carlo simulation and data fitting. The baseline concentrations of HbO and HbR were set to be 60uM and 40uM, respectively [28].

Blood flow calculation
For each frame of the laser speckle image, a spatial sliding window is used to calculate the contrast K as [Eq. (10)] [1]: where σ is the standard deviation and <I> is the average intensity in the sliding window. The relationship between contrast K and velocity v is [4]: where τ c is the decorrelation time, T is the exposure time, β is a factor related to detector size, speckle size and polarization.

rMRO 2 calculation
The rMRO 2 is calculated as [15]: where, v BF,0 represents the blood flow at the time 0 while ∆v BF represents the change in blood flow, here v BF was considered to be approximated as 1/ K 2 . γ R , γ T are vascular weighting constants which were set to be 1 [11].

Yeast-hemoglobin-intralipid phantom
The accuracy of hemoglobin concentrations and blood oxygen measurement was validated by the yeast-hemoglobin-intralipid phantom experiments. The phantom was a mixture of 0.3mL of chicken blood, 80mL phosphate buffer solution, 2mL of 20% intralipid and 0.8g yeast. The intralipid was used as scattering agent while the yeast was used to consume oxygen, and the initial concentration of HbT was approximately 5.8uM. The phantom was evenly mixed by a magnetic stirrer and pumped into a plastic tube placed on a piece of light-absorbing cloth by a peristaltic pump (BT100-2J, Longer, China). Firstly, the hemoglobin concentration of the phantom was changed by adding whole blood gradually. Chicken blood was added 15 times during the experiment and the increment was 0.3mL. The final volume of chicken blood was 4.8mL. The hemoglobin concentrations were calculated with our integral algorithm and the modified Lambert-Beer's law with a constant DPF corresponding to the initial concentration respectively.
In the second phantom experiment, the oxygen content of the phantom was changed by feeding excess air lasting for 6 minutes. The phantom was prepared in a similar manner as above, except that the chicken blood was increased to 3mL. Before the onset of the experiment, the evenly-mixed solution was stood alone for 10 minutes. The blood oxygen was calculated with two methods: the modified Lambert-Beer's law with the extinction coefficients at the peak wavelength of LEDs, and our integral algorithm shown in [Eq. (9)].

Arterial occlusion experiment
A six-minute arterial occlusion experiment was performed to observe hemodynamics changes during ischemia and reperfusion in human. A cuff was applied to the upper arm of an adult volunteer. After a two-minute rest, the pressure was increased to 220mmHg for two minutes, followed by a two-minute postrelease. Hemodynamic parameters of the hand were measured continuously with a sampling rate of 1 fps.

BFP and BVP observations
The experiment was carried out on four volunteers aged from 24 to 26. The volunteers placed their hands in the imaging area. Here the CCD1 was replaced with a monochrome camera with a high frame rate (1/2 inch, VGA, PCO Computer Optics, Germany). Each measurement lasted 40 seconds and the sampling rate was 50 fps. A region of interest (ROI) was chosen on the nail area of the right index finger. The blood flow and hemoglobin concentrations were averaged over the ROI and then analyzed in time series. A 0.8-6Hz Butterworth band pass filter was used to filter the signals.

Yeast-hemoglobin-intralipid phantom validation
The measured change of HbT concentrations when adding whole blood to the phantom was shown in Fig. 3. The standard deviation of the measured HbT concentrations at different ROIs of the phantom in a single experiment was calculated and plotted. The linearity of the results decreased with the increase of hemoglobin concentration when calculated with a constant DPF corresponding to the initial concentration. But, since the dependence of DPF on hemoglobin concentrations was considered, the measured hemoglobin concentrations calculated using our algorithm maintained a good linear response in the whole process. Figure 4 showed the change of hemoglobin concentrations as the blood oxygen changed. Figure 4(a) showed the results calculated using our integral method in [Eq. (9)] while Fig. 4(b) showed the results calculated with molar extinction coefficients at peak wavelengths of LEDs. Fig. 4(a) and Fig. 4(b) showed the same trend. During continuous air feeding, the concentration of HbO increased while the concentration of HbR decreased. After the air stopped, the concentration of HbO decreased and the concentration of HbR increased as the yeast consumed oxygen. When comparing the HbT concentration measured with the two methods, the results calculated with the integral algorithm were more stable, which meant more accurate results could be obtained using our algorithm, as shown in Fig. 4(c).

Arterial occlusion experiment
The spatiotemporal changes in rBF, ∆HbO, ∆HbR, ∆HbT, rMRO 2 at different times during the arterial occlusion were plotted in Fig. 5(a). The changes of rBF, ∆HbO, ∆HbR, ∆HbT over time in the ROI were shown in Fig. 5(b). And the change of rMRO 2 over time in the ROI was  illustrated in Fig 5(c). When the pressure was increased to 220mmHg, the blood flow was rapidly cut off, and the concentration of HbR continued to rise while the concentration of HbO continued to decrease. When the pressure was released, blood flow increased rapidly and then gradually fell back to the baseline. The concentrations of HbO and HbR quickly returned to the baseline, while the concentration of HbT raise by 10% and returned to the baseline. The rMRO 2 decreased quickly due to the decrease of blood flow perfusion at the beginning of the occlusion, then slightly increased because tissue oxygen delivery was lower than tissue oxygen consumption. Finally, the rMRO 2 increased rapidly because of the increasing perfusion after the postrelease.

BFP and BVP observations
The BFP and BVP observations of 4 subjects were shown in Fig. 6. The BFP signal and the BVP signal both showed heartbeats, as shown in Fig. 6(a). One period of the signals was taken out for analysis, as shown in Fig. 6(b). There was a time delay between BVP signals and BFP signals. The histogram distribution of the time delay in 40s was calculated for a single subject, which was concentrated at 40-100ms, as shown in Fig. 6(c). The individual average time delays of the four subjects were distributed between 57-137ms, while the total average time delay of the four subjects was about 90ms, as shown in Fig. 6(d).

Discussion
The relay lens system was composed of commercial achromatic lenses and meniscus lenses, which made the system simple in structure and low in cost. An infinite conjugate structure was adopted to reduce the beam aperture angle and make the relay distance more flexible. The maximum beam aperture angle in the relay unit was about 4.2°, which had little influence on the transmittance properties of optical filters. Therefore, the system is compatible with filter elements such as liquid crystal tunable filters (LCTF), which usually have stricter requirements for incident angle than normal filters have. By simply replacing the filters or increasing the number of the relay lens, the system can be extended to more applications, such as multispectral imaging and fluorescence imaging. Moreover, the relay system is also suitable for other imaging devices with C-Mount interfaces, such as endoscopes, colposcopes, and etc.
The relay distance of the system was mainly limited by the aperture of the relay lens units. As the relay distance increases, the image quality might decrease and vignetting might appear due to the aperture angle in object space of the relay system was smaller than that in image space of the camera lens. It is helpful to choose a short relay distance or relay lenses with large diameters to reduce image vignetting. Here a relay distance of 60 mm was selected considering the image quality and the size of the beam splitter.
The minimum F number of the imaging system was designed to be about 3, while the minimum working distance was designed to be 90 mm. Considering the Bayer array in the color CCD, if the acceptable diameter of confusion was thought to be 2 pixels (11 um), the minimum depth of field of the imaging system was about 8 mm and the depth of focus was approximately 0.9 mm. The depth of field can be increased by reducing the aperture or increasing the working distance.
To realize the stable and rapid acquisition of 50 fps, we chose a monochrome camera with PCI interface. The color camera was triggered by the monochrome camera. The effective image region was the overlapping field of view of the two CCD sensors. Image registration played an important role in this study. If the cameras are not moved or replaced, the registration can be performed only once to obtain the transformation matrix. For the convenience of the spatial image registration between two cameras, a monochrome and color version of the same USB3.0 camera could be used in the future.
Because of the aliasing spectrum response between blue channel and red channel as well as the wide spectrum distribution of the LED, there was a crosstalk between blue channel and red channel so that the received signals could be contaminated. The ratio of the two channels in 2.1 can represent the degree of crosstalk. Ideally, the red channel of CCD was expected to only response to the red LED, and the blue channel of CCD was expected to only response to the green LED, which means that the ratio of the two channels is expected to be infinite. The higher the ratio, the lower the crosstalk. This crosstalk can be reduced by choosing a camera with less aliasing spectrum response between channels.
Hemoglobin concentrations in different biological tissues vary widely. In traditional methods of estimating hemoglobin concentration, only the DPF corresponding to the baseline concentration was considered. In this case, the estimating results remained linear when tissue hemoglobin concentrations changes slightly. However, as the hemoglobin concentrations increase, the actual DPF will decrease due to the increase of absorption coefficient so that the linearity of measurement will be reduced. In this study, the dependence of DPF on HbT concentrations was considered in the calculation, so the measured HbT concentrations maintained better linearity.
During the yeast-hemoglobin-intralipid phantom validation experiment, although the actual HbT concentration remained constant when the phantom's blood oxygen changed, the measured HbT concentration increase slightly, as shown in Fig. 4. It was the typical crosstalk between blood oxygen and blood volume mentioned in previous studies [28]. Compared with the method using molar extinction coefficients at the peak wavelength of the LED, our integral method showed lower crosstalk. The accuracy was improved because our integral algorithm considered the integral effect of both the spectrum distribution of the LED and the spectrum responses of the camera. Reducing crosstalk between CCD channels as discussed above may help further improved the accuracy.
The imaging system was tested by the BVP and BFP observation experiment. A ROI on fingernail was selected because the nail bed contained an abundant network of capillaries, and the nail deck has good light transmittance. More stable signals can be acquired from the ROI on fingernails than on other regions covered with the complete epidermis structure. BVP signals reflected cardiac synchronous expansion and contraction of blood vessels with each heartbeat, while BFP signals represented the fluctuation of blood velocity during cardiac cycle [37,38]. Similar to the results in previous study, a time delay between BFP and BVP signals was observed. The time delay may be due to the fact that changes in blood flow velocity preceded changes in blood volume and have physiologic relevance for assessing microvascular flow and resistance [38,39].
The research had limitations. Biological tissues are not always flat. Accurate hemodynamic could hardly be obtained for tissues with curvatures exceed the depth of field of the system. Because an electrically tunable focusing lens was used in our imaging system, this can be further improved by using multi-focus image fusion [40].
In this study, only the blood volume was changed in the first phantom experiment while only blood oxygen was changed in the second phantom experiment. In Fig. 3, compared to the modified Beer-Lambert's law, the measured concentrations from the integral algorithm showed better linearity and larger value. But in Fig. 4, the measured concentrations from the both methods showed similar scale. It is probably because the effect of blood oxygen on the DPF was not considered in the integral algorithm. In the integral method, only the dependence of DPF on the HbT concentration was taken into account. However, oxygen content also affects DPF. Under the green LED illumination, HbO and HbR have similar absorption coefficients, so the DPF was mainly related to the HbT concentration. But the absorption coefficients of HbO and HbR are quite different under the red LED illumination, so the DPF depended on both HbT concentrations and oxygen content [32]. Two-dimensional fitting of DPF, blood volume and blood oxygen may help make further improvements.
In conclusion, a large field of view, multimodal imaging system was developed to measure multiple hemodynamic parameters. An infinite conjugate relay lens system compatible with the C-Mount camera lens was designed by using the commercial achromatic lenses and meniscus lenses, which extended the back focal length of the camera lens so that beam splitters and optical filters could be inserted into and multiple cameras could be connected. A electrically tunable focusing lens was introduced to allow rapidly focusing. An integral algorithm is proposed considering the dependence of DPF on hemoglobin concentrations and the integral effect of both the emission spectrum of light source and the spectrum response of detector. The system could be a helpful and extendable tool for exploring the macroscopic tissue hemodynamic responses.