1. Introduction

Cutaneous burn wounds are characterized by significant spatiotemporal shifts in tissue water content (TWC) both within and adjacent to the directly injured (contacted) region. Injury induced TWC variation has been characterized primarily with ex vivo wet/dry ratio techniques [1–3] and nuclear magnetic resonance (NMR) [4]. In vivo measurements of the edematous response to burn injuries have focused on tonometry, extremity dimeter (where applicable), and other morphology-based techniques that infer water content indirectly [5].

Histopathologic analysis of tissue excised from burn wounds in animal models have identified three, concentrically arranged (nominally), injury zones identified by the likelihood of self-healing. Proceeding from the injury center these are the zones of (1) coagulation, (2) stasis, and (3) hyperemia as shown in Fig. 1. In controlled animal models, the zone of hyperemia lies beyond the contact region, the zone of stasis typically lies adjacent to the contact region and can extend beyond the contact/no-contact border. Similarly, the zone of coagulation typically lies within the contact region and can extend beyond the contact/no contact border. Animal model studies suggest that the spatial extents of these zones are dynamic and correlate with wounds severity. Further, zone location and extent correlate with the distribution of abnormal (raised/reduced) TWC. Therefore, it is may be possible to characterize burn depth from quantification of wound edema and its spatiotemporal distribution.

 

Fig. 1 Burn zones organize as concentric like shells in the depth are hypothesized to appear as projections of the shells on the surface. Skin layers denoted A, B, C and damage zones denoted 1, 2, 3, 4. (a) Cross-sectional view: damage from superficial thickness wounds are limited primarily to the stratum corneum and epidermis, partial thickness burn wounds extend into the epidermis, and full thickness wounds extend into deep dermis. (b) Isometric view providing orientation between (a) and (c). (c) Top down view: The damage zones, and thus variations in surface TWC should present as regions concentric with the center.

Download Full Size | PDF

The edematous response of cutaneous tissue to thermal injury and its association with injury depth motivate the application of THz imaging and sensing as a diagnostics tool. The THz permittivity of water is large with respect to the all other non-aqueous tissue constituents [6]. Further, water is the only known tissue constituent that demonstrates significant dispersion at THz frequencies. These features allow THz systems to generate significant contrast from variations in TWC.

These concepts are demonstrated in Fig. 1. The depth dependent injury severity, as described by histology [7], is displayed in Fig. 1(a). Severity is organized as approximately concentric shells. The dotted lines indicate equal distance contours from the wound center. The surface projection of these concentric zones is displayed in Fig. 1(c) with Fig. 1(b) providing orientation. Injury severity is organized concentrically from the wound center. Consequently, tissue edema, and thus THz contrast should be organized about the burn wound center and along the tissue surface.

One major challenge in the study of THz imaging of cutaneous burn wounds are methods to relate THz contrast to visible contrast. The contact area, and thus injury center, are identified by burn specialist through visual appearance. Wound induction protocols often employ some form of contact and the contact area follows a controlled spatial extent (i.e. circular or rectilinear) with uniform heating within the perimeter. These protocols produce burn injuries with stark discoloration between contacted and non-contacted areas enabling rapid identification of wound landmarks and geometry.

Data in the literature suggest that the acute edematous response or cutaneous injury is large, rapid, and dynamic [1–3]. However, while the expectation is the edematous response is somewhat concentric to the wound center, the fact that the spatial distribution has not been characterized remains. This presents a problem to the analysis of THz imaging data since the majority of burn wound diagnostic imaging techniques are explored by referencing data to some geographical wound center and wound perimeter. The center and extent of the wound cannot be inferred directly from the THz imaging data thus the center must be characterized from a related data set and then mapped to the THz imaging data set.

This paper describes a methods and protocols by which the wound perimeter and center is located in visible images the burn wound. Further, the location of fiducial markers, arranged on the surface of the field-flattening window, and in the periphery of the imaging field are demarcated. These, markers generated significant contrast in both visible and THz spectra. Control points in the visible images were derived from manually drawn burn periphery contours. These points were used to locate the wound center in the THz image and spatially co-register all the THz and visible images. Automated segmentation was applied to the THz image data for unsupervised calibration of reflectivity values This set of operations produced registered THz data sets with a known, demarcated wound center and a shared, amplitude range. The processed image data was used to explore spatio-temporal variations in THz reflectivity in an animal burn injury model. Finally, the potential effects of dielectric window induced occlusion are discussed.

2. Image registration in THz medical imaging

Image registration is an essential task in diagnostic medical imaging as clinical decisions are often made using data sets from numerous, independent imaging modalities. This necessitates the identification of landmarks with contrast evident in data from all relevant modalities, and then creating a shared spatial domain where the mapped landmarks are collocated. A common baseline [1–3] for image registration is the visible appearance or visible light micrographs for the tissue of interest. Many clinical judgments, especially in the management of burn wounds, are based on visual inspection [8–10].

THz medical image registration is an active research topic. The visual appearance of diseased or injured tissue has been extensively characterized for the majority of applications under exploration. However, the THz imaging contrast within, and adjacent to, affected tissue is still poorly characterized in many applications. Sufficient analysis of spatiotemporal variation of THz contrast requires relating the observed spatial and temporal THz signatures to known visible image signatures. Historically, registration of THz images to visible images has been aided by optically clear, field flattening windows. These windows create a planar, specular surface amenable to THz data acquisition and establish a constrained, often circular field of view (FOV). If visible images of the tissue can be acquired without perturbation of window position, between THz and visible image acquisition, then approximate image registration can be achieved through manual shift, rotation, and dilation. Manual registration has also been enhanced through the addition of fiducial markers such as gold patterns on dielectric windows [11]. Ex vivo image registration is often aided by the extent of the tissue section not filling the THz imaging system FOV and thus the perimeter, defined by the large deviation in THz contrast between sample and no sample can be mapped to the pronounced sample border apparent in the visible images [12].

There are limited demonstrations in the literature regarding registration between visible and THz image sets of burn wounds. Spatial correlations between visible and THz imagery is apparent in experiments from Tewari et al. [13] and Bajwa et al. [14] but no image registration was performed.

This work details methods that (1) registered companion THz and visible image data of cutaneous burns in rat models, (2) registered all image from all rats to the same shared space, and (3) calibrated THz data using pixel population derived from automated, threshold-based segmentation. This shared space is called wound space and defines pixels by their distance from the wound center. Following THz images co-registration, the spatiotemporal dynamics of burn wound THz reflectivity in tissue, defined by its position relative to the wound center, was analyzed.

3. System

A visible image and optical component ray diagram of the reflective THz imaging system are displayed in Fig. 2(a) and Fig. 2(b), respectively. The THz source, receiver electronics, and system design have been previously reported [15, 16] and are summarized below.

 

Fig. 2 (a) Block diagram of the reflective THz imaging system with a mode locked laser (MLL) driven photoconductive switch (PCS) and Schottky diode detector. (b) THz optical system consisting of four off-axis parabolic (OAP) mirrors. The windowed burn injuries were placed coincident to the target plane defined by the overlap of OAPs 2 and 3. (c) Visible image of the THz system in the animal surgery room. An SLR camera was used to acquire concomitant visible pictures of the optical field of view. Infrared (IR) data was also collected but not used in this study.

Download Full Size | PDF

A GaAs photoconductive switch (PCS) was pumped by a frequency doubled, mode locked laser (MLL) with the following parameters: τ_p = 100 fs, f_rep = 46 MHz λ₀ = 780 nm, P_ave = 120 mW. The PCS free space output was collimated and focused by off-axis parabolic mirrors (OAP 1, 2 in Fig. 2(b)). The rat abdomen was placed coincident with the target plane with the beam incident at 14°. The reflected THz radiation was collimated and focused by a second set of OAPs (OAP 3, 4 Fig. 2(b)) to the feedhorn of a WR1.5 waveguide mounted Schottky diode detector. The rectified THz pulses were amplified (38 dB) and coupled to a gated receiver driven with a reference RF pulse time locked to the THz pulses. Pixels were generated with a 1 ms integration time yielding a peak of SNR of 40 dB. THz Images were generated by raster scanning the rat model beneath the fixed, focused THz beam using x- and y-axis stepper motors. Acquisition time for a 60 mm x 60 mm field of view (FOV) with 0.5 mm x 0.5 mm pixels (120 x 120 = 14,400 pixels) was ~7 minutes.

The PCS power spectral density coupled with the schottky diode spectral responsivity yielded an operational band centered at 525 GHz with ~125 GHz of bandwidth. The spectral densities combined with the optics enabled a diffraction limited spot size of ~1 mm at a 36 mm standoff distance.

An image of the THz imaging system installed in the animal operating room is displayed in Fig. 2(c). The rat was placed on the indicated “animal mount”. The THz imaging system, digital SLR, and midwave infrared (not reported in this study) imaging camera were arranged in close proximity enabling the enabling precise, repeatable, stepper motor translation between the imaging modalities.

4. Experimental protocol

The in vivo study was approved by the UCLA Institutional Review Board (#2009-094-03b). 10 male Sprague Dawley (SD) rats weighing 500- 700 g (Harlan laboratories, Hayward, CA) were used as pre-clinical models. The total number of animals (n = 8) was determined by statistical requirements of a paired, summed rank test (Wilcoxon t-test) . The computation assumed a power of 0.9, a significance level of 0.05, a simulated mean THz reflectivity signal difference of 0.5% for each scan, and a computed difference standard deviation of 0.424%. Preliminary data was utilized from previous studies [13, 15, 16] for the power calculations. Two extra rats (total of 10 subjects) were included in the study to account for random death. One animal in the partial thickness burn wound group died leaving four animals in the partial thickness group and five in the full thickness group.

Each rat was anesthetized using isoflurane (4% and 1% for induction and maintenance, respectively), shaved to expose a 5 cm x 5 cm area of bare abdominal skin, and aseptically cleaned with three alternating scrubs of Betadine and 70% alcohol. Rats were then placed in supine position on a mount with internal heating to maintain normal body temperature during imaging. Each rat was subcutaneously injected with Buprenorphine (0.05 mg/kg) 1 hour prior to imaging.

An imaging window was lowered onto the shaved abdomen of the animal to flatten the imaging field and minimize image contrast due to the irregular surface geometry. The window was constructed using 50.8 mm outer diameter, 38.1 mm inner diameter, 3 mm thick polished brass ring and a suspended 12.7 thick µm Mylar film (Fig. 3(a)). Additionally, triangular sections of blue tape were adhered to the brass ring to provide control points for image registration (Fig. 3(a)). Visible images of the uninjured abdomen were captured with a SLR camera and reference scans of the uninjured abdomen covering a 60 mm x 60 mm field of view (FOV) were acquired with the THz imaging system. The window was removed and a 5 mm x 19 mm rectangular brass brand was applied to the abdomen with light, constant contact pressure for 10 seconds for both burn types. Partial thickness burn wounds were induced at 130°C and full thickness burns at 200°C. The thermal parameters were verified with histology in a control population prior to the start of this study. Finally, as a preventative measure against dehydration, all rats were intraperitoneally injected with a high dose of pharmaceutical grade Lactate Ringer’s Solution (60ml/kg) following burn induction. A visible image of a partial thickness burn beneath the Mylar window is shown in Fig. 3(a). The burn is outlined by the red dotted contour.

 

Fig. 3 Experimental protocol. (a) Visible image of a partial thickness burn under the Mylar window suspended by the brass ring. The white arrows indicate the blue triangular tape sections the burn injury perimeter is denoted by the red dotted contour. (b) Excised, burned abdomen tissue with the metal tissue resection guide. The burn perimeter is outlined in red and the section locations are indicated by the black dotted lines. Tissue sections spanned both uninjured and injured tissue to aid histological assessment. (c), (d). Representative histology of sections from a partial thickness and full thickness burn respectively Reference axes in the bottom left corner or the image align the x-axis with the long dimension of the burn injury.

Download Full Size | PDF

Following burn induction, the brass ring mounted Mylar window was repositioned on the burned abdomen and concomitant visible and THz images were continuously acquired at the following time points (in minutes): [15, 30, 45, 60, 75, 105, 135, 165, 195, 225, 255, 285, 315, 345, 375, 405, 435]. Note that the window was engaged throughout the 435 minutes of scans. A visible image showing a partial thickness burn wound and brass mounted Mylar window are displayed in Fig. 3(a). Following the scans, the rats were revived and returned to the vivarium. Trimethoprim-Sulfamethoxazole (antibiotic) was administered orally to prevent possible infection after burn injury. Individual follow-up with visible and THz images were acquired at 24 hr, 48 hr and 72 hr post burn induction (Not included in this work, see [17]).

Animals were euthanized after acquisition of the 72 hr scan. Three sections of skin along the rostral-caudal axis were harvested from the left, center and right regions of the burn using a section guide (Fig. 3(b)) and transferred to 10% formalin solution. The tissue sections spanned both uninjured and injured tissue, thus providing a control area to compare against the burned area. The skin samples were then submitted for histological embedding, sectioning, and staining with hematoxylin and eosin (H&E) (Fig. 3(c)-(d)). Light microscopy was used to examine to observe structural tissue damage to assign burn severity subject to standard criteria [10].

5. Image registration and wound space

Image registration was performed by rigid affine transformations consisting of translation and rotation (Eq. 1) driven by manually selected control points represented by P₁, P₂, P₃, and P₄ in Fig. 4. Shearing was not necessary due to accurate control of rat position relative to the THz and visible light imaging systems. Scaling was performed separately though image resampling. The inner diameter of the brass ring (38.1 mm) was used as a known distance standard and bicubic interpolation was used to down-sample the SLR images by a factor of ~8 and up-sample the THz images by a factor of 4 resulting in equal sampling density grids of 62.5 µm X 62.5 µm.

 

Fig. 4 Visual description of the two-step registration for visible and THz imaging data. (a) Raw camera frame (e.g. of rat 2) showing a burn with manually defined perimeter contour and manually selected control points P₃ and P₄ at the apex of the blue tape. (b) Best fit of a 14 mm line segment within the manually segmented contour. The endpoints of the line are demarcated control points P₁ and P₂. (c) Registration of Rat 2 to Rat 1 by collocating points P₃ and P₄ from both images. (d) Control point displacement from the raw camera frame to the wound space frame. (e) Raw THz image frame (e.g. of rat 2) with manually selected control points P₃ and P₄ at the apex of the blue tape. (f) THz frame registered to the wound space through colocation of P₃ and P₄. (g) Superposition of wound space center on registered THz frame. The insets i₁-i₄ show the displacement between P_1- P₄ and P^'₁- P^'₄.

Download Full Size | PDF

Pixels representing the intersection of the blue tape triangles and brass window inner circumference were demarcated by hand, in both the visible and THz images of the burn wounds. These control points are denoted P₃ and P₄ in Fig. 4(a) and Fig. 4(d). Next, a closed contour was hand drawn on the visible images of the burn wound (Fig. 4(b)). A 14 mm line segment (brand length – brand width = 19 mm – 5 mm), was fit within the hand drawn burn wound perimeter curve such that all the distance vectors between the line segment and perimeter contour were minimized (Fig. 4(c)). This line represented the burn wound center and thus the 0-distance contour. The endpoints of the line segment, denoted P₁ and P₂ in Fig. 4 were demarcated as control points.

One of the nine rat models (herein referred to as rat 1) was randomly chosen as the reference frame to which the data from the other rats (herein referred to as rats 2-9) were registered. First, intramodality, interpatient registration was performed with Eq. (1) to map the 0-distance contours in the visible images of rats 2-9 to the 0-distance contour in rat 1. P₁ defined the translation (t_x, t_y in Eq. (1)) and the angle between the (P₁,P₂) line segment of rat 1 and the (P₁,P₂) line segment of rats 2 – 9 defined the rotation (θ in Eq. (1)) about P₁. This process is described in Fig. 4(c) and Fig. 4(g) where the mapped control points are denoted by the prime symbol: P₁ → P^'₁ P₂ → P^'₂, P₃ → P^'₃, P^'₄ → P₄.

Affine transformations were also applied for intermodality, intrapatient registration where the THz images were registered to the transformed visible image sets. The mapped locations P^'₃ and P^'₄ were demarcated as the baseline. The displacement between visible image P^'₃ and THz image P₃ defined the translation (t_x, t_y) and the angle between line segment (P^'₃, P^'₄) and line segment (P₃, P₄) defined the rotation (θ).

This process resulted in the center (0- distance contour) of every acquired THz and visible image lying coincident in shared space termed “wound space”. Wound space enabled grouping and analysis of pixel values with respect their distance from the brand contact center and hence wound center. Note that 17 THz images and 17 visible images (digital SLR) were acquired for each rat in the 435 minutes following injury. While significant temporal variation in THz image contrast was observed throughout this time, no changes were discernable in the visible images. Therefore, the visible image acquired at 15 minutes was used to extract the wound center that was mapped to each set of 17 THz images.

6. Image calibration

THz pixel amplitude calibration was accomplished by segmenting the THz images with amplitude thresholding techniques and then using the populations of segmented pixels to compute reference values. The circular window mount underfilled rectilinear FOV and ensured that each THz image contained corner regions corresponding to zero-signal which were used as the “0%-reflectivity” reference. The flat brass ring that surrounded the abdomen served as a “100%-reflectivity” reference and enabled line-by-line calibration of system drift as well as detection of imaging plane tilt due to either physical tilt, or temporal system drift.

6.1 Segmentation

Air: (1) A first-pass estimate of air pixels was obtained via threshold-based segmentation. (2) A mask consisting of pixels near the corners was constructed using pixels closest to the corners (closest 40%) that also lie within a region based on the corners of a circle inscribed by a square. (3) The intersection of the threshold-based segmentation and the corner-based segmentation was computed and tabulated as “air pixels”. (4) The median of the “air segmentation” was used as the 0%- reflectivity for each THz image. The segmented air regions are demarcated ‘iii’ in Fig. 5(a) and the pixel population is displayed in Fig. 5(b).

 

Fig. 5 Image segmentation and calibration. (a) Cropped wound space THz frame and pixel masks denoting (i) brass, (ii) abmomen, and (iii) air. (b) Pixel histograms associated with regions i, ii, and iii. Note the increased histogram spread in the brass ring pixels due to tilt and shot noise. (c) Evaluation of least squares plane correction on a flat reference reflector. The left column demonstrates the algorithm performance on a leveled surface. The right column is the same reflector tilted 7 degrees about the y-axis as defined in Fig. 2(b).

Download Full Size | PDF

Brass: (1) A first-pass guess of brass pixels was determined via threshold-based segmentation. (2) The connected components of the segmentation were found and the largest components were added until 30% of the area was accounted for. (3) Similar to the air segmentation, a mask consisting of areas close to the corners of the image was generated with the additional constraint that the brass regions cannot intersect the corners. (4) Lastly, the segmentation is morphologically eroded to prevent the edges of the brass from contributing to the calibration routine. The segmented brass region is demarcated “metal pixels” and the median of this population was used as 100%-reflectivity for each THz image. The segmented brass region is demarcated ‘ii’ in Fig. 5(a) and the pixel population is displayed in Fig. 5(b). The increased histogram spread with respect to the segmented air pixels was primarily due to tilt. Shot noise also likely contributed (see section 9).

Abdomen: Abdominal segmentation was found by finding the largest circular mask having minimum overlap with the brass segmentation. (1) The (x,y) coordinates of the brass ring mask center were computed. (2) For every point within 1 cm of the central point, a circle of fixed radius was constructed (chosen to be slightly smaller than the inner radius of the brass ring), and the overlap between this circle and the brass ring was computed. (3) (x,y) coordinates that produced the minimum overlap were used to construct the largest possible circle that did not overlap the brass ring. The segmented abdomen regions are demarcated (ii) in Fig. 5(a) and the pixel population is displayed in Fig. 5(b).

6.2 Least squares plane correction

Following segmentation, a least squares plane in R³ was fit to the pixels within the segmented brass region of each THz image. $I_0\left[x_k,y_k\right]$ in Eq. (2) is the amplitude of the kth pixel in the brass ROI located at row $y_k$ and column $x_k$ and v₁, v₂, and v₃ are constants extracted from the least squares fitting routine.

The fitted plane was subtracted from the image to give the segmented brass a mean value of ~1 (Eq. (3)) and then image contrast was linearly stretched such that the mean value of the brass remains 1 and the mean air reflectivity is mapped to 0 (Eq. (4)).

6.3 Segmentation and calibration evaluation

The center of a polished brass plate was positioned coincident with the virtual rotation point of a two-axis goniometer. The goniometer assembly was leveled and mounted such that the virtual point was coincident to the THz imaging system focal point. THz images of the target at angular settings ranging from −7 to 7 degrees, in steps of one degree about the x and y axes defined in Fig. 2(b), were imaged.

A THz image of the characterization target at 0 degrees tilt (referred to as the “baseline” image) is displayed in the top left panel of Fig. 5(c). The 7 deg. tilt about the y-axis created the largest deviation across the segmented brass and is displayed in the top right panel of Fig. 5(c). Contours were drawn using the techniques described in Section 6.1 and the least squares plane correction algorithm, described in Section 6.2, was applied. The corrected 7 degree tilt image is displayed in the right panel, in the middle row of Fig. 5(f). The perturbation due to applying correction to an already level data set (0 deg. image) is displayed in the left panel of Fig. 5(c). Both demonstrate good uniformity within the segmentation contour.

The performance of the least squares algorithm was evaluated by computing the per-pixel mean square error (MSE) and total average image MSE between the algorithm corrected images in the row labeled “tilted” in Fig. 5(c) and the baseline image, within the demarcated ROI. The MSE maps are displayed in the bottom row of Fig. 5(f) where each image is paired with its own color bar/map. The plane correction algorithm applied to the baseline image imparts minimal change in image intensity variation with an average MSE of 2.5x10⁻⁷. The distribution of MSE reveals a small bias towards the bottom left corner of the target, suggesting there may be a very minor tilt in the baseline target not controllable with the leveling technique. The seven-degree tilt corrected images yielded an average MSE of 8.2x10⁻⁴ with peak MSE at the edges of the ROI. This error is nearly an order of magnitude smaller than the 0.54% reflectivity differential used in the statistical power calculation and indicates sufficient performance.

7. Registered and calibrated THz imaging results

7.1 Image display

The results of image registration and calibration are demonstrated in Fig. 6 for data from a full thickness burn acquired at 435 minutes. The registered visible and THz images are displayed in Fig. 6(a) and (b) respectively. The shifted and rotated original FOV is outlined in white. The arrows in both figures denote the wound center contours. The image rotation and shift operations insert 0’s outside of the original FOV as indicated by the figure text.

 

Fig. 6 Image superposition. The registered, calibrated THz image (b) is superimposed on the registered visitble image (a) to create the figure in (c). The calibrated colormap is denoted in (d) and the transparency applied to this colormap is denoted in (e) with the transparency percentage listed in parenthesis. The vertical line in (e) demonstrates amplitude varying transparency.

Download Full Size | PDF

The THz image was overlaid on the visible image with a “hot” colormap and a transparency determined by the reflectivity range of pixels subtended by the abdomen mask. Maximum abdomen reflectivity of was set to white with a transparency of 0 while minimum reflectivity was set to black with a transparency of 1. Linear mappings were used to compute the reflectivity and transparency inside these ranges. Figure 6(c) shows the results of overlaying Fig. 6(a) with Fig. 6(b) where the visible image is obscured by THz data in the high THz reflectivity areas. The pairing of reflectivity and transparency are demonstrated in the colorbars in Fig. 6(d), (e).

7.1 THz images of partial and full thickness burns on Day 0

Representative THz time-series images from a partial thickness burn wound and full thickness burn wound versus time are displayed in Fig. 7 (a) – (j) and Fig. 7(k) – (t), respectively. Each burn series consists of two rows. The top row depicts THz contrast in a hot color map where black-to-white corresponds to minimum and maximum abdomen tissue reflectivities. The bottom row superimposes the THz reflectivity maps onto the registered visible light image of the burn wound where the colormap was modified to a black-to-cyan pallet with the same transparency mapping demonstrated in Fig. 6.

 

Fig. 7 Burn wound time series imaging results. (a) – (e) THz images of a partial thickness burn. (f) – (j) Partial thickness THz images superimposed on the registered visible frame. (k) – (o) THz images of a full thickness burn. (p) – (t) Full thickness THz images superimposed on the registered visible frame. THz contrast is distinct for each burn wound severity. In time-series THz images of a partial thickness burn, the contact area shows a drop in TWC and an edematous front superior to the burn over time. In contrast to the partial thickness wound the contact area of a full thickness burn does not display a significant drop in TWC. Additionally, the contact area is surrounded by a ring of TWC which runs concentric and coincident with the burn contact zone.

Download Full Size | PDF

The first column of each set contains the control images for both rat models. The unburned abdominal tissue appears fairly uniform in the visible and THz images with natural variation consistent to what was observed in previous studies. The next column contains the first images acquired following thermal insult. The time point is labeled 15 minutes; this is the result of the ~7 minutes required to acquire the image acquisition and ~8 minutes to move the rat model from the branding station to the imaging system, reposition the window, and reconnect the gas anesthesia. Scans were performed every 15 minutes up to 75 minutes and then every 30 minutes, terminating at 435 minutes. During this time the rat models remained under the imaging system with the Mylar window engaged.

The partial thickness burn wound demonstrated an initial, THz reflectivity increase throughout the FOV with a localized decrease in THz reflectivity approximately spanning the contact area. The appearance is similar at 105 minutes although the reflectivity increase, immediately superior to the contact region has transitioned into an high reflectivity front that runs along the entire top edge of the contacted area. Conversely, the contact area maintains a general drop in reflectivity around the periphery and a thin area of elevated reflectivity; comparable to pre-burn values. 195 minutes shows continuous upward movement of the high reflectivity front superior to the wound and the emergence of a high reflectivity front inferior to the wound. By 435 minutes, the edematous fronts have reached the top and bottom of the FOV. A small localized increase in THz reflectivity organized just inferior to the contact area and through the contact area center.

The full thickness burn wound demonstrated an initial significant THz reflectivity increase in the periphery and an overall drop in THz reflectivity in the contact area. At 105 minutes, the peripheral THz reflectivity begins to while the contact THz reflectivity begins to increase. This behavior is continued through the 195-minute mark. By 435 minutes, distinct regions are visible in and around the contact region. In contrast to the partial thickness wound, the contact area displays a significant increase in THz reflectivity. Additionally, the contact area is surrounded by a ring of decreased THz reflectivity which runs concentric and coincident with the contact zone.

8. Analysis of dynamics

The variation in THz reflectivity as function of location with respect to (1) the wound center and (2) time elapsed following injury was computed using a 2-point, pixel-wise derivative. The derivatives were then averaged to produce the maps in Fig. 8(a). This process is described by Eq. (5) where I_k[m,n] is THz image k acquired at time t_k, [m, n] are row and column location of a particular pixel, N = 17 is the total number of images acquired in the post burn, 435 minute trial, and $\overline{\dot{I}\left[m,n\right]}$ is the average derivative map. Figure 8(a) superimposes nine aggregate abdomen masks which are each the logical AND operation of the 17 abdomen masks per rat. The global, shared overlap is indicated by the shaded area and the co-registered 0-distance contour is represented by the black contour. The average derivative maps were cropped to the shared overlap mask detailed in Fig. 8(a) thus limiting the analysis to pixel locations present in all 17 x 9 = 153 THz images.

 

Fig. 8 THz reflectivity dynamics during the 435 experiment duration. (a) Perimeters of all eroded abdomen masks from all 9 rats. The shared overlap (logical AND operation of all masks) is highlighted by the grey shaded area. (b) – (e) pixel-wise average rate of change of THz reflectivity in %/min for the partial thickness burn injuries. (f) – (j) pixel-wise average rate of change of THz reflectivity in %/min for the full thickness burn injuries. The minimum average rate of change occurs within the contact area for 8 of the 9 animals.

Download Full Size | PDF

The partial thickness burn spatio-temporal maps are shown in Fig. 8(b) – (e) and the full thickness burn wounds are in Fig. 8(f) – (j). Each image has its own colormap in units of % change in THz reflectivity per minute. The displayed colorspace is specific to each image and stretches across the full range of the average computed derivatives.

The spatio-temporal variation maps illustrate three interesting phenomenon. (1) The contact areas for 8 of the 9 rats, (b)-(e) and (g)-(j), demonstrated significantly reduced changes in THz reflectivity compared to the periphery. The minimum average reflectivity change in the 0-distance contour over 7 + hours is 0.03% while the maximum in the periphery varies from 0.15% to 1.01% corresponding to increase of 3x to 34x. Further, the minimum computed derivative occurs within the contact area for the 8 indicated rat models.

(2) The transition between low variation in the contact area and high variation in the periphery is generally rapid (~0.3%/mm) and the border between these areas appears to span parts of the zone of stasis for the majority of the wounds and contrast between the low.

(3) There is a visible rostral-caudal bias in the spatiotemporal variation maps. The abdomen region superior to the 0-distance contour (closer to the head) demonstrated an average, space integrated variation of ~0.13%/min while the region below the contour was a factor of ~3 down at ~0.04%/min. The contact region demonstrated an average time integrated variation of ~0.02%/min.

The spatiotemporal variation of full thickness burn wound in Fig. 8(f) is distinctly different than its cohort. The approximate zone of stasis featured some of the lowest observed changes in THz reflectivity within that animal’s data while the zone of coagulation featured increased rates of change over the periphery. Nothing remarkable was noted in histology and the animal did not significantly differ from the others in terms of age, weight, etc. The source of the discrepancy is still uncharacterized.

It’s important to note that the very first derivative was computed with time t₁ = 15 minutes and time t₂ = 30 minutes. Since the perimeter of the burn was used to inform image registration, it was not possible to use the pre-burn image and thus, the change in reflectivity between uninjured and the t = 15 minutes could not be characterized. Therefore, one can conclude that majority of THz contrast variation, correlated with the edematous response, was observed in tissue peripheral to the contact region, in the period spanning 15 minutes to 435 minutes. It is possible that the initial, acute response in the first 15 minutes following injury demonstrated changes in THz in the contact region that exceeded changes in the subsequent time.

Overall, we believe that the observed spatiotemporal variations support the potential utility of regions in the periphery of the contact area. Further, these observations suggest the edematous response in areas adjacent to and/or outside the visibly injured region may hold information relevant to the status of the burn wound.

9. Dielectric window engagement and shot noise

Continuous dielectric pressure/contact is a potential confounder and recent studies have reported on the resulting effects of tissue occlusion in healthy human subjects [18]. The effects of window-induced occlusion on THz imaging contrast was explored by acquiring THz images of uninjured abdominal skin, in two anesthetized rats, with continuous window engagement over a 7hr period. System drift was calibrated using the methods described in Section 6. Figure 9(a)-(d) contains representative images at time points 15 min., 90 min, 240 min, and 420 min respectively from one of the rat models. The data is expressed with a “hot” color palate that maps the minimum THz reflectivity (air) to black and the maximum THz reflectivity (brass) to white. The images series demonstrate a general decrease in reflectivity throughout the abdomen. Note that the start time was defined by window placement, not burn induction hence the difference in acquisition times points compared to Fig. 7

 

Fig. 9 Exploration of occlusion effects from application of dielectric window. (a)-(d) Representative THz images of uninjured tissue beneath the Mylar window from four different time points. The hot color map transforms black to the global minimum THz reflectivity and white to the global maximum. Reference axes are in the bottom left corner of each panel. (e) Ribbon plot of THz reflectivity across the field of view corresponding to the location of the superimposed cyan contours in (a) – (d). The circular markers in (e) map to the circular markers in (a) – (d). (f) Per-pixel standard deviation map of the field of view demonstrating signal variance over the 7 hour study.

Download Full Size | PDF

The contours superimposed on Fig. 9(a)-(d) were positioned to span abdomen areas of the highest and lowest observed variation. The time evolution of pixels along these contours is displayed as a ribbon plot in Fig. 9(e). The x-axis in Fig. 9(e) correspond the x-axes in Fig. 9(a)-(d) and the circular markers aligned with t-axis in Fig. 9(e) correspond to the circular maker endpoints in Fig. 9(a)-(d). The ribbon plot shows an initial drop in the reflectivity of the periphery of ~5% hour in the first 30 minutes and an additional drop of ~3% over the subsequent 390 minutes (indicated by the black arrow). The central region (0 mm – 18 mm) of the abdomen maintains a reflectivity within ~1.5% of its initial value with negligible changes observed over the 8 hour experiment.

The overall variation was summarized by constructing a time series (19 time points) at every pixel location in the FOV and computing the standard deviation of each time series. This standard deviation map is displayed in Fig. 9(f). The periphery demonstrated a reflectivity variation of ~6% over 8 hours while the center of the abdomen demonstrates a variation of less than 1%. Also notable is brass region signal variation which ranges from 1% - 2% and often exceeds that observed in the abdomen center. This behavior is consistent with shot noise; an effect observed in the brass reflectivity PDF of Fig. 5(c).

The global drop in reflectivity over time is consistent with [18] although the difference systems, protocols, and animal vs human models precludes comparison of specific, computed reflectivity values.

9.1 Discussion

The 12.7 μm thick, n = 1.6 Mylar window is optically thin at THz frequencies was selected because it created contrast nearly identical to that expected from a non-contact, windowless setup. However, due to Mylar’s flexibility, it is possible that the pressure exerted by the brass ring mount is greater towards the edge of the FOV as compared to the center due to the relative flexibility of the 12.7 μm film compared to the inner diameter edge of the brass ring. The signal variation on the periphery is higher than the center although it monotonically decreases throughout the experiment duration; in kind with the FOV center. The potential for window contact based confounders cannot be ruled out but. An occlusion trial on rats, constructed similar to [18], is warranted.

Brass was selected as the reference reflector because it resists oxidation, is easy to clean/polish, and is relatively dispersion free in the THz band [19]. When combined with air, brass regions provide a convenient two-point reflectivity calibration. However, the entirety of the brass region displayed higher temporal variance than the lowest temporal variance areas of the abdomen. Further, the air regions displayed the lowest temporal variation. The positive correlation between signal amplitude and signal variance strongly supports the presence of shot noise. Replacing the brass ring with a rigid, well characterized dielectric will produce less noisy calibration and improve the results of the calibration routines.

10. Conclusions

An experimental protocol and associated image processing techniques were presented that enabled spatial registration of visible and THz images of cutaneous burn wounds. The methods relied on manual selection of fiducial marker derived control points in both the visible and THz image sets as well as manual segmentation of the burn wound perimeter in the visible image. The segmentation was reduced to a set of control points representing the center of the contact area and, thus, the center of the wound space. Affine transformation were used to place visible and THz images of cutaneous burn wounds in a shared space such that the wound centers were all collocated. The major benefit of this technique is that THz images can be registered without relying on THz frequency wound contrast which cannot be known a priori. This construct allows evaluation of wound contrast as a function of distance from the contact area across a large population of animals. Automated segmentation was applied to the THz image data for unsupervised calibration of reflectivity values.

THz images of partial and full thickness burn wounds in nine animal models were acquired over 7.25 hours and registered to a shared wound space. Contrast variation maps were computed for each animal by computing derivatives with respect to time for every pixel location and image acquisition time point. The absolute values of the derivatives were averaged to create a picture of tissue contrast dynamics. The maps generally showed THz tissue reflectivity variations were greater on the rostral side of the wound as compared to the caudal side. Additionally, the contact area tissue displayed significantly less change in THz reflectivity than the periphery tissue.

An imaging trial of uninjured abdomen tissue over 7 hours revealed a monotonic drop in THz reflectivity as a function of time suggesting that occlusion, due to window contact, may be a confounder.

The results elucidated significant spatiotemporal variations in injured tissue THz reflectivity. The periphery, non-contacted tissue, exhibited significant variation in contrast and supports the exploration of diagnostic utility of tissue beyond the visibly affected region. This finding is unique and may set THz apart from competing biophotonic modalities. The presented protocol would benefit from reduced contact pressure, faster acquisition time, and redesigned fiducial makers to enable better capture of THz reflectivity variations immediately following thermal insult and reduce the effect of occlusion on the acquired THz data.