Patent Application
20200340798
The presently disclosed subject matter relates generally to imaging. Particularly, the presently disclosed subject matter relates to systems and methods for imaging based on multiple cross-sectional images acquired at different angles to enhance resolution.
Optical coherence tomography (OCT) is a cross-sectional, micron-scale imaging modality based on coherence gating for depth resolution that has become the clinical standard of care for pathological diagnosis and treatment monitoring in several medical specialties. The axial resolution is determined by the coherence length of the light source and can be sub-micrometer. However, in a tradeoff necessitated by diffraction, conventional implementations of OCT typically accept an inferior lateral resolution, on the order of 10 μm or greater, in order to obtain a long depth of focus for bulk tissue imaging. Previous studies have addressed this tradeoff using techniques such as digital refocusing and beam-shaping, but in all cases the lateral resolution was still limited by diffraction.
Another limitation of conventional OCT is that images are distorted due to spatially inhomogeneous refractive index (RI) distributions in tissue. Previous attempts to correct such RI-induced distortions required a priori information on both the sample geometry and RI values, which severely limits the generalizability of such techniques.
In view of the foregoing, there is a need for improved OCT systems and techniques that limit image distortions and provides improved image resolution and contrast.
Having thus described the presently disclosed subject matter in general terms, reference will now be made to the accompanying Drawings, which are not necessarily drawn to scale, and wherein:
The presently disclosed subject matter provides systems and methods for imaging based on multiple cross-sectional images acquired at different angles. According to an aspect, a method includes acquiring multiple cross-sectional images of an object at different angles. The method also includes registering the acquired cross-sectional images. Further, the method includes reconstructing an enhanced resolution image of the object based on the registered images. The method also includes displaying an image of the object based on the enhanced resolution image.
According to another aspect, an OCRT system includes a conventional imaging system, to generate 2D or 3D cross-sectional images of a sample. The system also includes a mechanism acquire these images at angles up to 360 degrees. Furthermore, the system includes an algorithm that registers the 2D/3D images, reconstructs a higher resolution 2D or 3D image of the sample, and generates a 2D or 3D RI map of the sample. Further, the system includes a computing device that implements the aforementioned algorithm and displays the result.
The following detailed description is made with reference to the figures. Exemplary embodiments are described to illustrate the disclosure, not to limit its scope, which is defined by the claims. Those of ordinary skill in the art will recognize a number of equivalent variations in the description that follows.
Articles “a” and “an” are used herein to refer to one or to more than one (i.e. at least one) of the grammatical object of the article. By way of example, “an element” means at least one element and can include more than one element.
“About” is used to provide flexibility to a numerical endpoint by providing that a given value may be “slightly above” or “slightly below” the endpoint without affecting the desired result.
The use herein of the terms “including,” “comprising,” or “having,” and variations thereof is meant to encompass the elements listed thereafter and equivalents thereof as well as additional elements. Embodiments recited as “including,” “comprising,” or “having” certain elements are also contemplated as “consisting essentially of” and “consisting” of those certain elements.
Recitation of ranges of values herein are merely intended to serve as a shorthand method of referring individually to each separate value falling within the range, unless otherwise indicated herein, and each separate value is incorporated into the specification as if it were individually recited herein. For example, if a range is stated as between 1%-50%, it is intended that values such as between 2%-40%, 10%-30%, or 1%-3%, etc. are expressly enumerated in this specification. These are only examples of what is specifically intended, and all possible combinations of numerical values between and including the lowest value and the highest value enumerated are to be considered to be expressly stated in this disclosure.
Unless otherwise defined, all technical terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this disclosure belongs.
As referred to herein, the term “computing device” should be broadly construed. It can include any type of device including hardware, software, firmware, the like, and combinations thereof. A computing device may include one or more processors and memory or other suitable non-transitory, computer readable storage medium having computer readable program code for implementing methods in accordance with embodiments of the present disclosure. A computing device may be, for example, retail equipment such as POS equipment. In another example, a computing device may be a server or other computer located within a retail environment and communicatively connected to other computing devices (e.g., POS equipment or computers) for managing accounting, purchase transactions, and other processes within the retail environment. In another example, a computing device may be a mobile computing device such as, for example, but not limited to, a smart phone, a cell phone, a pager, a personal digital assistant (PDA), a mobile computer with a smart phone client, or the like. In another example, a computing device may be any type of wearable computer, such as a computer with a head-mounted display (HMD), or a smart watch or some other wearable smart device. Some of the computer sensing may be part of the fabric of the clothes the user is wearing. A computing device can also include any type of conventional computer, for example, a laptop computer or a tablet computer. A typical mobile computing device is a wireless data access-enabled device (e.g., an iPHONE® smart phone, a BLACKBERRY® smart phone, a NEXUS ONE™ smart phone, an iPAD® device, smart watch, or the like) that is capable of sending and receiving data in a wireless manner using protocols like the Internet Protocol, or IP, and the wireless application protocol, or WAP. This allows users to access information via wireless devices, such as smart watches, smart phones, mobile phones, pagers, two-way radios, communicators, and the like. Wireless data access is supported by many wireless networks, including, but not limited to, Bluetooth, Near Field Communication, CDPD, CDMA, GSM, PDC, PHS, TDMA, FLEX, ReFLEX, iDEN, TETRA, DECT, DataTAC, Mobitex, EDGE and other 2G, 3G, 4G, 5G, and LTE technologies, and it operates with many handheld device operating systems, such as PalmOS, EPOC, Windows CE, FLEXOS, OS/9, JavaOS, iOS and Android. Typically, these devices use graphical displays and can access the Internet (or other communications network) on so-called mini- or micro-browsers, which are web browsers with small file sizes that can accommodate the reduced memory constraints of wireless networks. In a representative embodiment, the mobile device is a cellular telephone or smart phone or smart watch that operates over GPRS (General Packet Radio Services), which is a data technology for GSM networks or operates over Near Field Communication e.g. Bluetooth. In addition to a conventional voice communication, a given mobile device can communicate with another such device via many different types of message transfer techniques, including Bluetooth, Near Field Communication, SMS (short message service), enhanced SMS (EMS), multi-media message (MMS), email WAP, paging, or other known or later-developed wireless data formats. Although many of the examples provided herein are implemented on smart phones, the examples may similarly be implemented on any suitable computing device, such as a computer.
As referred to herein, the term “user interface” is generally a system by which users interact with a computing device. A user interface can include an input for allowing users to manipulate a computing device, and can include an output for allowing the computing device to present information and/or data, indicate the effects of the user's manipulation, etc. An example of a user interface on a computing device includes a graphical user interface (GUI) that allows users to interact with programs or applications in more ways than typing. A GUI typically can offer display objects, and visual indicators, as opposed to text-based interfaces, typed command labels or text navigation to represent information and actions available to a user. For example, a user interface can be a display window or display object, which is selectable by a user of a computing device for interaction. The display object can be displayed on a display screen of a computing device and can be selected by and interacted with by a user using the user interface.
Embodiments of systems and method disclosed herein are referred to “optical coherence refraction tomography” or “OCRT,” which can use multiple OCT cross-sectional images (“B-scans”), which can be 2D or 3D, acquired at a diversity of angles (in 2D or 3D) to reconstruct isotropic, high-resolution, cross-sectional images with the superior axial coherence gating of conventional OCT extended to the lateral dimension. Alternatively, if the OCT system has a higher lateral resolution than axial resolution (such as in spectroscopic OCT (SOCT)), the superior lateral resolution may be used to improve the axial resolution up to the lateral resolution. Taken to the other extreme, if the OCT system is completely defocused or unfocused such that there is effectively no lateral resolution, OCRT still theoretically increases the lateral resolution up to the axial resolution limit. Here, “OCT” and “conventional OCT” are to be interpreted broadly to include all of its variants, such as, but not limited to, time-domain OCT (TD-OCT), spectral-domain OCT (SD-OCT), swept-source OCT (SS-OCT), full-field OCT (FF-OCT), line-field OCT (LF-OCT), optical coherence microscopy (OCM), spectroscopic OCT (SOCT), pump-probe OCT (PP-OCT), microscope-integrated OCT (MI-OCT), and stimulated Raman scattering OCT(SRS-OCT), all of which can have the shortcomings of anisotropic spatial resolution and refraction-induced spatial distortions. More generally, instead of OCT one may also use any imaging system with anisotropic spatial resolution, such as ultrasonic imaging, photoacoustic imaging, and their variants. Even more generally, instead of OCT, we may also use any imaging modality that suffers from sample-induced image distortions. However, herein the term “OCT” or “conventional OCT” is used. OCRT systems and methods disclosed herein can be used to address these shortcomings of conventional OCT. In the absence of refraction, light rays travel in straight lines and the reconstruction procedure is analogous to that of X-ray computed tomography (CT). However, due to refraction and optical path delays in the sample, OCT images acquired at multiple angles are distorted differently and hence need to be dewarped and registered. Thus, in correcting for sample-induced refraction, OCRT is unique for estimating the spatially resolved RI distribution of the sample, which is aligned with the high-resolution reconstruction. Other multiangle RI tomography techniques for thick samples accounted only for path delay (not changes in ray direction), required access to the other side of the sample, and were not applied to biological samples. Focal shift-based OCT methods for depth-resolved RI measurement typically require experimental repositioning of the focus and physical depth scanning in the sample or reference arm. However, OCRT can be implemented using, for example, a suitable OCT system by rotating the sample or by use of a beam scanning protocol or a combination thereof in accordance with embodiments disclosed herein. The multiangle measurements need not be acquired sequentially, but can be acquired in a multiplexed fashion, such as by having multiple OCT systems viewing the sample at different angles, or by having a modified OCT system that simultaneously collections multiple angular views. Such strategies may be combined with sample rotation, beam scanning, or both.
The system 100 may include an OCT scanning device 106 configured to scan an object 108, such as biological tissue, and to output image data representative of the scanned object 108. The image processor 100 may receive the output image data and store the data in memory 102. A user interface 110 (e.g., a keyboard and/or mouse) may be used by an operator (e.g., ophthalmologist, radiologist, etc.) to enter input. The image processor 100 may process the image data in accordance with embodiments disclosed herein. The processed data may be stored in memory 102 and may be displayed as an image and/or text on a display 112 of the user interface 110.
As OCRT can be described as a Fourier synthesis technique from multiangle illumination, it is useful to compare it to CT. A single 1D projection in CT can be considered an image resulting from the convolution of the 2D scene with a point-spread function (PSF) that is infinitely wide in projection dimension (x) and narrow in z. For example,
In some embodiments, OCRT and CT may be considered conceptually similar (See e.g.,
Now referring to
While RIs of biological tissue in the X-ray regime are close to 1, in the optical regime a spatially varying RI distribution distorts OCT images by changing the path lengths and the directions of the rays, preventing rigid-body registration among the multiangle B-scans. If the RI distribution are known, the rays may be propagated accordingly to dewarp the image prior to FBP application to generate a high-resolution reconstruction. To infer the unknown RI distribution, in OCRT the inverse problem may be solved using, for example, the ray equation as the forward model, which in 2D (x and z) is
where nA(x,z) is the RI distribution parameterized by A and s is the position along the 1D ray trajectory. The parameterization can be a sum of a regularly spaced grid of Gaussian kernels such that nA(x,z) is differentiable everywhere and minimizes the effects of the “staircase” artifacts stemming from discretization onto a Cartesian grid. In particular, an index distribution given by the Nadaraya-Watson kernel parameterization is one example choice, which avoids having to use finite differences to compute spatial gradients, as analytic expressions are available. Other suitable parameterizations of the RI distribution may be used. Additionally, other methods for solution of the forward ray propagation problem (such as other analytical or computational approaches which are known in optics) may be used. Furthermore, in addition to methods for ray propagation, methods for wave propagation and photon diffusion may be used.
Because in general there is no closed solution to the ray equation, numerical integration (fourth-order Runge-Kutta) was employed. To specify optical path delay, a constant step size may be used for the numerical solver, which can be scaled by the inverse of the RI value at the current position of the ray. To initiate ray propagation, the initial conditions are to be specified, which correspond to the initial ray positions and directions. For example, assuming telecentric scanning and uniform sampling of A-scans within each B-scan, the initial directions are all parallel, and the initial positions are equally spaced for each angle. In this way, the A-scans can be propagated from each B-scan from each angle as individual rays, which cause the overall B-scan image to dewarp. In theory, if the optimal RI distribution is found, the images would be perfectly registered. It can be assumed that the group and phase indices are approximately the same in biological tissue (minimal chromatic dispersion), such that a common index distribution governs both path trajectory and delay. Thus, herein the terms group (refractive) index and RI are used interchangeably. However, it is possible for OCRT to reconstruct the group and phase refractive indices separately by treating them as separate parameters in the optimization algorithm.
To provide feedback on the accuracy of nA(x,z) to aid its optimization, a differentiable metric may be used that quantifies the degree of joint registration among all the B-scans. Here, an intensity-based metric is described, by which the mean squared error (MSE) between the raw B-scans and a forward prediction of the B-scans may be computed based on the estimated high-resolution reconstruction. However, other metrics may be used, such as those used by point-set-based registration techniques. However, other metrics may also be used. For example,
A modified version of gradient descent called Adam may be used to minimize this regularized MSE, which jointly registered the B-scans and estimated both the RI distribution and the high-resolution, isotropic image. Other gradient-based optimization methods may be used such as stochastic gradient descent, Adagrad, Adadelta, and Nesterov momentum, among many others. The gradient of the intensity-based registration metric may be computed through the numerical differential equation solver with respect to forward model parameters using TENSORFLOW™, a software library that employs automatic differentiation, a widely used technique in the deep learning community requiring only the specification of the differentiable forward model from which the gradient is computed through recursive application of the chain rule. Other deep learning software may be used, such as PYTORCH, CHAINER, DYNET, among many others. More generally, any framework that can implement gradient-based optimization may be used. Although gradient-based optimization is described here, gradient-free optimization algorithms such as evolutionary algorithms and other stochastic optimization methods may be used. In summary, the OCRT optimization procedure simultaneously registers all B-scans, generates an undistorted, isotropic, high-resolution reconstruction, and a co-aligned estimate of the RI distribution of the sample.
To validate the isotropic resolution of OCRT, 560-nm polystyrene beads embedded in a 2% (w/v) agarose gel inside a glass microcapillary tube were imaged, as shown in
To quantify the resolution of OCRT, the beads of both the reconstruction and the averaged volume were localized. Also, a 2D Gaussian function with axes oriented with the Cartesian coordinate system was fit, where z and x are the axial and in-plane lateral dimensions, respectively (for the reconstruction, the x and z axes are arbitrary). The resulting distributions of the PSF width fits in x and z (See
To validate the RI distribution estimates, another bead phantom was made substituting the agarose gel for polydimethylsiloxane (PDMS).
To demonstrate the performance of our method, OCRT was imaged and applied to a variety of ex vivo biological samples, including several mouse organs (vas deferens, femoral artery, bladder, and trachea), a human donor cornea, and a marsh crane fly (Tipula oleracea) leg. All samples were inserted into microcapillary tubes for convenience during imaging, and hematoxylin and eosin (H&E)-stained histological sections of neighboring tissue samples were obtained for comparison (also Masson's trichrome stain was obtained for the artery sample). Overall, it was observed that lateral resolution was significantly improved in the OCRT reconstruction across all samples imaged, and that the features in the RI maps matched features in the OCRT reconstructions.
In the first vas deferens sample (panels “a”-“h” of
In the mouse femoral artery sample (see panels “q”-“x” of
In the mouse bladder sample of panel “a”-“h” of
Finally, the OCRT reconstructions of the tibia and femur cross sections of a crane fly (j) exhibit superior resolution of the outer walls (cuticle) and the ultrafine features that surround them. The cuticle, apparent as a thin orange structure in the H&E, manifests at the bottom of femur as two thin dark lines in the OCRT reconstruction. It is noted that although the larger hair-like structures (setae) present in the B-scan and OCRT reconstruction are not visible in the H&E, likely lost during tissue processing, the much finer setae lining the wall in the H&E are still visible in the OCRT reconstruction. The circumferential orientation of the internal longitudinal muscle fiber is also much more apparent in the OCRT reconstruction than in the B-scan. Furthermore, the RI map indicates that the lumen of the tibia was filled with air (n≈1) while the lumen of the femur was filled with water. Although this is disparity is likely due to the sample preparation procedure, it is a feature not at first apparent from the B-scan. However, upon closer inspection, we note that the region indicated by the green arrow in
As disclosed herein, OCRT is provided and may use OCT B-scans from multiple angles to simultaneously reconstruct an undistorted, isotropic, high-resolution image and estimate the RI distribution. In particular, the superior axial resolution of OCT is extended to the inferior lateral resolution to form an isotropic PSF. Resolution enhancement and RI estimates of OCRT are demonstrated not only in sub-resolution scattering bead phantoms, but more importantly in several biological samples.
As an example, OCRT is advantageous over other techniques that address the depth-of-focus tradeoff because it can replace the lateral resolution dependence on diffraction with the axial resolution dependence on coherence. For example, recent work in applying OCT to the extreme ultraviolet regime has yielded axial resolutions on the order of 10s of nm, but with lateral resolutions three orders of magnitude lower. Since OCRT does not require high-NA objectives, much longer working distances associated with low NAs are possible, which make in vivo imaging more practical. Moreover, OCRT can tolerate lateral aberrations, in particular those induced by the sample. While the resolution of OCRT can be degraded by chromatic dispersion induced by the sample, in principle they can be corrected computationally. Another advantage of OCRT over other inverse optimization approaches is that it does not require phase stability, which for near infrared wavelengths typical of OCT requires nanometer scale stability of both the scanner and the sample over the acquisition time. On the other hand, OCRT, which employs intensity images, only requires stability on the order of the axial resolution, and even this constraint may in principle be relaxed through computational corrections during the registration step. Another advantage is the speckle reduction due to angular compounding of independent speckle patterns.
In accordance with embodiments, OCRT imaging systems and methods disclosed herein provide a general framework that leverages multiangle OCT images to computationally synthesize both an enhanced-resolution reconstruction and an RI map. Substantial improvements of OCRT over OCT in various biological samples are described, and it was found that there were structures readily apparent in OCRT reconstructions but missing in conventional OCT images.
OCT images shown in the figures were acquired using a commercial spectral domain OCT system (Bioptigen Envisu R4110)MR SDOIS) at a 20-kHz A-scan rate, with an incident power of 1 mW, an 820-nm center wavelength, and a nominal axial resolution of 1.2 μm in air and lateral resolution of 8.5 μm. The B-scans consisted of 500 A-scans over a field of view of 1.5 mm; each A-scan had 2048 pixels with a maximum imaging depth of 2.22 mm in air. This is just one example of an imaging system that acquires distorted, lower resolution images at multiple angles; however, as mentioned earlier, other imaging systems may be used (e.g., other variants of OCT, ultrasound imaging, photoacoustic imaging, etc.).
In one use case, for each sample, 20-averaged OCT B-scans or volumes every 6° across 360° were acquired. For the polystyrene bead phantoms only, out-of-plane averaging of B-scans was performed across 18 frames spaced 3 μm apart to account for slight tube misalignment, due to which the beads close to the wall of the tube away from the axis of rotation, because of their small size, were not present in every B-scan. These y-averaged B-scans were used as the raw data for OCRT reconstruction. For all biological samples this was not an issue and so the B-scans were not averaged across the y direction.
All samples were inserted into capillary microtubes (provided by Drummond Scientific Microcaps) with an inner diameter of 797.6 μm and immersed in phosphate buffered saline (PBS). 560-nm polystyrene beads (provided by Thermo Scientific) were embedded in 2% agarose gel (provided by Sigma Aldrich) or PDMS (provided by Dow Corning). Tube insertion is an example technique to facilitate multiangle image acquisition, in this case through sample rotation.
The animal organs used in this study were from mice (C57BL/6 wild-type). The organs were fixed in 10% neutral buffered formalin (NBF) and kept at 4° C. for 24 hours. After fixation, all organs were micro-dissected using an upright Zeiss dissection scope. The dissected tissues were inserted into the glass microtubes using dissection forceps. In some instances, the tissues were placed into the glass tube following aspiration with a 23G truncated needle on the opposite side of the tube. This procedure allowed the specimen to gently slide into the inner wall of the glass microtubule. To avoid dehydration after fixation, all specimens were transferred with the glass tube in PBS. The remaining organ tissue was used for paraffin block preparation followed by H&E staining or trichrome staining for histological examination. Histological images were acquired at 10× or 20× magnification using an Olympus microscope and digitally white balanced.
N=8 tubes were sampled and their inner diameters estimated using OCT by measuring the difference between the luminal reflections, and obtained 797.3±0.9 μm (standard error of the mean), in very good agreement with the manufacturer's specification of 797.6 μm. The inner diameter pathlength of the bead samples (with either agarose or PDMS embedding media) were estimated using the same method for all angles and volume (y) slices, across which the standard deviation values in
To acquire images at multiple angles, a custom-built, inverted rotation stage (Thorlabs) was designed and made to mount a single tube vertically. The tube had two orientation degrees of freedom and two translational degrees of freedom perpendicular to the axis of rotation, which we use to manually align the tube. The tube was immersed in water or PBS in a cuvette, which remained stationary during sample rotation. The entire above setup was mounted on xyz-translation stage as shown in
Optimization was conducted on TensorFlow 1.8 using Python 2.7 on an Intel Core i7-3930K with 48 GB of RAM. Because each iteration required 10s of GB of memory, the CPU version of TensorFlow was used. For all samples we ran gradient descent for 200-500 iterations, on the order of a 2-3 minutes per iteration. Filter optimization was run for an additional 100 iterations. The same hyperparameters were used for all biological samples.
Here, an example image dewarping model based on the ray equation is provided, although other suitable image deformation models may be used. A 2D ray equation is given by Equation 1, but because s is the position along the 1D ray trajectory, we may use one of the two equations, and use the relationship
Thus, we used a single nonlinear, second-order differential equation:
As discussed herein, because there is in general no closed solution to the ray equation, a fourth-order Runge-Kutta (RK4) numerical integration was used on Equation 2 to obtain the ray trajectory, x(z). To specify changes in optical path length, the RK4 step size can be divided by the RI at the current position. In principle, a ray may be propagated for each A-scan, but in practice an approximation may be used and propagated as a smaller number of rays. Because the RK4 solver is an iterative chain of element wise tensor multiplication operations, it is differentiable.
Equation 2 requires the spatial gradient of nA at any arbitrary position (x,z). Furthermore, because gradient descent may be used for optimization. Second-order, mixed partial derivatives of nA may be available (i.e.,
for all α∈A). One such parameterization is the Nadaraya-Watson parameterization using the Gaussian kernel, given by
which is a weighted sum of Gaussian kernels of uniform width σ and heights given by matrix A, which in our experiments had dimensions 128 by 128 covering a field of view of ˜1.6 by 1.6 mm. The kernel positions may be given by xi, zi, which were regularly spaced on a 2D grid in experiments. An analytical expression for the gradient of Equation 3 at an arbitrary position is also available. This parameterization may be chosen because sharp index boundaries were approximated as smooth transitions, especially if the RK4 step size is larger than the resolution of the index distribution. It was found that propagating rays at a boundary between two uniform media with this approximation agreed with Snell's law. Another reason was that we encountered “staircase” discretization artifacts, whereby an interface oblique to the coordinate axes would be locally parallel or perpendicular if some form of higher order interpolation is not used.
To ease the computational burden, several approximations may be made. For example, it was found that setting the RK4 step size to be the axial length of one pixel (in air) of an A-scan was too computationally demanding. Instead, the following may be made: let the step size cover p=8 pixels (˜8 μm) and use linear interpolation made to fill in the intermediate positions within the step. Likewise, it was found that propagating one ray per A-scan (˜380) was also too computationally demanding, so the rays may be downsampled to a fixed k=35 rays. The coordinates of the intermediate rays may be spaced evenly between neighboring propagated coordinates.
Since Equation 3 in practice may be too computationally demanding, when the index and index gradient at (xi,z1) are sought, only a neighborhood of 7 by 7 kernels may contribute to the weighted sum (rather than every single kernel). As mentioned, it is noted that the first two approximations may affect the spatial resolution of the RI map.
To aid in optimization of the RI distribution, before A is allowed to vary arbitrarily, optimization with respect to parameters describing the tube geometry may be implemented. For example, there may be six parameters: the inner and outer radius of the tube, the medium index, the tube wall index, and the xz coordinates of the center of the tube. The inner and outer radii may be independently calibrated using OCT, and estimated the center of the tube in a separate optimization. Thus, only two RI values could be optimized. This may allow for a better initial guess for A. This procedure may be applied for a calibration sample, such as the polystyrene bead sample embedded in agarose gel, and values similar used to the resulting parameters may be used to seed the optimization for all biological samples.
The A-scans are backprojected along those trajectories according to the RI map rather than along straight lines as part of the image registration. Because fewer rays (k=35) than A-scans (˜380) might be propagated and each RK4 step size may cover p=8 axial pixels (Sec. 1.3), for all the other points the two closest rays may be identified to be propagated (one to the left, one to the right of the points) and on those rays the two closest depths (one above, one below) that were multiples of p (i.e., a depths at which the RK4 solver stepped). These four corner points describe the smallest rectangle (whose vertices were propagated points) that surrounded the points, which were evenly spaced with respect to these four points using a bilinear interpolation scheme (see
To do the backprojection, a tensor R of dimensions nx2×nxz×ny may be defined, choosing nxz such that the pixel size may be chosen appropriately to match the expected resolution of the reconstruction. In experiments, a single averaged B-scan was obtained from each angular view, so ny was set to 1. The values corresponding to the individual spatial coordinates for the B-scan pixel positions (which have been rotated, warped, and interpolated based on the propagated rays) may be distributed to fixed neighborhoods of pixels, weighted by the Gaussian kernel distance of each pixel to the spatial coordinate being projected. Other interpolation schemes may be used to reduce pixel-level artifacts.
Because the rays are no longer necessarily parallel after image dewarping, some parts of the final image could have more dense contribution (e.g., due to ray focusing); to account for this, normalization may be implemented by element-wise division by the coverage of the backprojection, which can be computed by backprojecting tensors of is whose sizes match those of the B-scans. In this way, the values accumulated in pixels in the reconstruction via backprojection that were visited by rays more times would be divided by larger numbers. This normalization scheme may be applied first for each B-scan, and finally for the full reconstruction.
Since backprojecting the raw B-scans may result in strong emphasis on the low spatial frequencies (
The intensity-based registration metric may be the mean squared error between the B-scan data and the B-scan prediction generated by the forward model from an iterative estimate of the high-resolution reconstruction.
Given the current estimate of the RI distribution, the rays may be propagated, and the B-scan pixel values were interpolated and backprojected as described herein. Let R denote the current estimate of the high-resolution reconstruction result. To generate B-scan predictions, the coordinates of the same propagated rays were used for backprojection to gather the reconstruction intensity values. In particular, for each coordinate along the rays, the intensities of a fixed neighborhood of pixels may be gathered and their mean computed, weighted by the Gaussian kernel distance to the pixels.
To make this B-scan prediction more accurate, optimizable parameters may be included that describe for example: 1) the angular dependence of backscattering, 2) the anisotropic point-spread function (PSF) of OCT, and 3) the axial attenuation (in this order).
OCT analyzes backscattered light whose intensity varies as a function of angle between the illumination and the orientation of the interface. Namely, if we image a flat surface, the backscattered intensity at normal incidence can be expected to be the highest and at near parallel incidence to be the lowest. Moreover, the angular dependence of the backscattered intensity depends strongly on the roughness of the surface (e.g., the backscatter intensity may be independent of angle for a sufficiently rough surface).
Since the angular dependence of any specific sample is unknown, it may be treated as an optimizable parameter by discretizing the angular range of −90° to +90° into for example d=30 levels, so there were d parameters to optimize, which may be used to rescale the B-scan intensity values pixel-wise. To estimate these d parameters, the reconstruction estimate R may be convolved with d Gabor filters oriented at those d angles and for each pixel in the resulting convolution the index corresponding to the angle with the largest Gabor response may be found, resulting in an image of local orientations (one of d indices). Next, the same ray trajectories may be used to gather (using rounded coordinates) the orientations, resulting in an image of the same size as the B-scan of orientation indices, each pixel of which was assigned a value according to the orientation index (one of d). Finally, the B-scan prediction may be multiplied by this image. Example optimized angular backscatter profiles for biological samples are given in
After the B-scan prediction is modified according to local orientation, it is convolved with a 1D Gaussian kernel, such that blurring only occurs in the lateral dimension.
Finally, a model may be incorporated to account for A-scan attenuation due to backscatter, absorption, and system sensitivity falloff. Accounting for attenuation may be necessary to prevent the RI from changing excessively from trying to account for registering A-scans at depths where there was no longer any signal. There was no attempt to separate the contributions to the axial attenuation, but rather nonparametrically modeled the total attenuation; as such, the parameters optimized for this purpose were not readily interpretable, but other parameterizations may be used to isolate the relative contributions. The following are defined:
where α, β, and γ are optimizable parameters. Bpsf is used in the exponential because light incident obliquely on a surface might result in lower backscatter, and thus the backscatter signal could underestimate the attenuation due, for example, to a specular reflection that missed the collection optics. The y exponent allows this attenuation model freedom to deviate from Beer-Lambert's law to account for other sources of attenuation (e.g., falloff). More generally, the exponential may be any parameterizable function of Bpsf, such as a Taylor expansion.
The MSE is then computed for each B-scan across all angles, which was then minimized iteratively.
Individual xz spatial translations for each B-scan may be allowed, to account for residual B-scan misalignment not explainable by RI. However, it was found that the amount of freedom allowed for the translations may affect the estimated RI values; in particular, it is possible for the optimization to move the B-scans away from the center and stretch the B-scans to compensate. Thus, regularization may be imposed on these translations as described herein.
Masked versions of L2 regularization and total variation (TV) regularization may be used on A. Instead of letting every pixel contribute equally to the regularization terms, the contributions may be weighted using a mask. In general, the outer boundary can in principle be determined because the shape of the first boundary is not distorted by refraction. Hence, L2 regularization may be used to constrain the region outside of the object being imaged to be close to the index of the immersion medium (e.g., the group index of water at 820 nm is ˜1.342). Similarly, the TV regularization may be restricted to the inside of the object (e.g., to preserve boundary between the immersion medium and the sample).
Two regularization terms may be introduced for the xz translations (δxθ, δzθ), which are individually defined with respect to the coordinate system of each B-scan. One is a simple sum of squares deviation:
Σθδxθ2+δzθ2 (5)
The other counteracts expansion away from the center by penalizing B-scans separated by 180° for having similar z shifts:
Σθ(δzθ+δzθ+180°)2 (6)
This term encourages the relative distance between the two opposing B-scans to be close to 0 and hence counteracts the reconstruction from expanding or shrinking.
Once all the differentiable terms contributing to the loss are defined, which is a scalar, TensorFlow may be used to compute the gradient of this loss with respect to all the aforementioned optimizable parameters. Subsequently, this loss may be iteratively minimized using the Adam optimizer.
The ray gathering step uses only a small neighborhood of pixels and thus if the B-scans are grossly misregistered the optimization may get stuck in local minima. To combat this, a multiresolution approach may be used. For example, decreasing nxz of R increases the spatial coverage of each pixel, effectively increasing the local search range. Thus, the process may start a small nxz reconstruction size and periodically increase it during optimization. A multiresolution approach was used for the human cornea samples, starting with nxz that was ⅛ of the final reconstruction size and switching to the final size halfway through optimization.
After the aforementioned parameters (RI map, spatial shifts, attenuation, etc.) are optimized, all parameters may be frozen that may alter the ray trajectories and iterative optimization may be performed of the backprojection filter using the same intensity-based metric. This allowed the reconstruction to be more consistent with the B-scans by reducing potential filter-induced artifacts. For example, it was found that using the Ram-Lak or the square root of the Ram-Lak filter introduced a halo due to low energy at low frequencies. Although in principle it is possible to optimize the filters jointly with the registration parameters, it was empirically found that using a filter that attenuates higher frequencies (e.g., square root of Ram-Lak) can promote registration (likely in a similar way that a multiresolution resolution approach benefits from starting with low resolutions).
In cases where the axis of rotation of the sample do not perfectly correspond to the center of the sample, each B-scan may be spatially translated so that the initial boundaries were registered. Fortunately, the initial boundary, in this case the tube outer wall boundary, is often easy to segment because it was the first bright reflection in the images. Then, let Q be the collection of all points corresponding to the initial boundary from all B-scans, let (x0, z0) be the coordinates of the virtual center of rotation that we wish to optimize such that the initial boundaries become registered, and let L be the axial-lateral aspect ratio. Assuming the absolute distance specified for the axial dimension of A-scans is trusted, L may be allowed to rescale the x coordinates. Next, without making any assumptions about the shape of the initial boundary, the full 360° was divided into angular slices {Δθj}j=0,1,2, . . . ′ where Δθj denotes the jth angular slice. The idea is that if the initial boundaries from multiple B-scans are well registered, the distances of initial boundary points from all B-scans to the center that reside in Δθj should be very close to each other, and thus we used the variance of those distances as a metric:
L
initial(x0,z0,L)=Σj var {∥((x−xc)L,z)−(x0,z0)|2|(x,z)∈(Δθj∪Q)}, (7)
where xc is a constant whose only purpose is to aid the optimization procedure by decoupling L and x0.
OCRT may require careful calibration of the xz aspect ratio of the B-scans, because the final reconstruction formed from rotated B-scans no longer has this distinction. It was found that the estimated absolute RI values were sensitive to this calibration, and that one effective way to calibrate the aspect ratio was to use the axial dimension of OCT to measure the diameter of the tube and then during registration of the circular initial boundary constrain the diameter of the tube to the measured diameter. That is, we can add an extra term to Equation 7,
λ(mean {|((x−xc)L,z)−(x0,z0)|2|(x,z)∈Q}−rmeas)2, (8)
where rmeas is the measured radius of the tube and λ is a coefficient tuning this term's importance. This method was particularly effective because the same “ruler” (i.e., A-scan depth) was used both to measure the diameter and to define the axial coordinates of the B-scans, such that even if the A-scan depth was off by a scale factor, it would cancel out. This of course depended on the tube being perfectly circular, which we found to be a reasonable assumption from the tube maintaining the same distorted shape in the B-scans upon rotation. For more general samples (i.e., not embedded in a tube), precise methods would be needed to calibrate this aspect ratio.
In addition to the human donor cornea sample presented in panels “a”-“h” of
In
In experiments, an additional mouse bladder sample were also imaged, which is in a relaxed state (the bladder sample in the main text was distended). Lateral resolution enhancement is clear in the OCRT reconstruction, which increased the visibility of the smooth muscle layers and lamina propria, consistent with the histology.
The present subject matter may be a system, a method, and/or a computer program product. The computer program product may include a computer readable storage medium (or media) having computer readable program instructions thereon for causing a processor to carry out aspects of the present subject matter.
The computer readable storage medium can be a tangible device that can retain and store instructions for use by an instruction execution device. The computer readable storage medium may be, for example, but is not limited to, an electronic storage device, a magnetic storage device, an optical storage device, an electromagnetic storage device, a semiconductor storage device, or any suitable combination of the foregoing. A non-exhaustive list of more specific examples of the computer readable storage medium includes the following: a portable computer diskette, a hard disk, a RAM, a ROM, an erasable programmable read-only memory (EPROM or Flash memory), a static random access memory (SRAM), a portable compact disc read-only memory (CD-ROM), a digital versatile disk (DVD), a memory stick, a floppy disk, a mechanically encoded device such as punch-cards or raised structures in a groove having instructions recorded thereon, and any suitable combination of the foregoing. A computer readable storage medium, as used herein, is not to be construed as being transitory signals per se, such as radio waves or other freely propagating electromagnetic waves, electromagnetic waves propagating through a waveguide or other transmission media (e.g., light pulses passing through a fiber-optic cable), or electrical signals transmitted through a wire.
Computer readable program instructions described herein can be downloaded to respective computing/processing devices from a computer readable storage medium or to an external computer or external storage device via a network, for example, the Internet, a local area network, a wide area network and/or a wireless network, or Near Field Communication. The network may comprise copper transmission cables, optical transmission fibers, wireless transmission, routers, firewalls, switches, gateway computers and/or edge servers. A network adapter card or network interface in each computing/processing device receives computer readable program instructions from the network and forwards the computer readable program instructions for storage in a computer readable storage medium within the respective computing/processing device.
Computer readable program instructions for carrying out operations of the present subject matter may be assembler instructions, instruction-set-architecture (ISA) instructions, machine instructions, machine dependent instructions, microcode, firmware instructions, state-setting data, or either source code or object code written in any combination of one or more programming languages, including an object oriented programming language such as Java, Smalltalk, C++, Javascript or the like, and conventional procedural programming languages, such as the “C” programming language or similar programming languages. The computer readable program instructions may execute entirely on the user's computer, partly on the user's computer, as a stand-alone software package, partly on the user's computer and partly on a remote computer or entirely on the remote computer or server. In the latter scenario, the remote computer may be connected to the user's computer through any type of network, including a local area network (LAN) or a wide area network (WAN), or the connection may be made to an external computer (for example, through the Internet using an Internet Service Provider). In some embodiments, electronic circuitry including, for example, programmable logic circuitry, field-programmable gate arrays (FPGA), or programmable logic arrays (PLA) may execute the computer readable program instructions by utilizing state information of the computer readable program instructions to personalize the electronic circuitry, in order to perform aspects of the present subject matter.
Aspects of the present subject matter are described herein with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the subject matter. It will be understood that each block of the flowchart illustrations and/or block diagrams, and combinations of blocks in the flowchart illustrations and/or block diagrams, can be implemented by computer readable program instructions.
These computer readable program instructions may be provided to a processor of a computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks. These computer readable program instructions may also be stored in a computer readable storage medium that can direct a computer, a programmable data processing apparatus, and/or other devices to function in a particular manner, such that the computer readable storage medium having instructions stored therein comprises an article of manufacture including instructions which implement aspects of the function/act specified in the flowchart and/or block diagram block or blocks.
The computer readable program instructions may also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable apparatus or other device to produce a computer implemented process, such that the instructions which execute on the computer, other programmable apparatus, or other device implement the functions/acts specified in the flowchart and/or block diagram block or blocks.
The flowchart and block diagrams in the Figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present subject matter. In this regard, each block in the flowchart or block diagrams may represent a module, segment, or portion of instructions, which comprises one or more executable instructions for implementing the specified logical function(s). In some alternative implementations, the functions noted in the block may occur out of the order noted in the figures. For example, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowchart illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardware-based systems that perform the specified functions or acts or carry out combinations of special purpose hardware and computer instructions.
While the embodiments have been described in connection with the various embodiments of the various figures, it is to be understood that other similar embodiments may be used, or modifications and additions may be made to the described embodiment for performing the same function without deviating therefrom. Therefore, the disclosed embodiments should not be limited to any single embodiment, but rather should be construed in breadth and scope in accordance with the appended claims.
This application claims priority to U.S. Patent Application No. 62/623,044, filed Jan. 29, 2018, and titled SYSTEMS AND METHODS FOR OPTICAL COHERENCE REFRACTION TOMOGRAPHY, the content of which is incorporated herein by reference in its entirety.
This invention was made with government support under grant number DGF-1106401 awarded by the National Science Foundation (NSF). The government has certain rights to this invention.
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/US19/15540 | 1/29/2019 | WO | 00 |
| Number | Date | Country | |
|---|---|---|---|
| 62623044 | Jan 2018 | US |
