WO2021174200A1 - Systems and methods for composing high dimensional microscopy images in subpixel resolution - Google Patents

Abstract

The present disclosure describes systems and methods for composing high dimensional fluorescence microscopy images with subpixel resolution. Several embodiments of the invention are directed to analyzing a plurality of images taken over multiple dimensions (e.g., X-, Y-, and Z-dimensions). Additional embodiments are directed to composing high dimensional fluorescence microscopy images taken across multiple cycles, where each cycle represents a series of images that are tiled in each of the X-, Y-, and Z-dimensions.

Description

SYSTEMS AND METHODS FOR COMPOSING HIGH DIMENSIONAL MICROSCOPY

IMAGES IN SUBPIXEL RESOLUTION

CROSS-REFERENCE TO RELATED APPLICATIONS

[0001] This application claims the benefit of U.S. Provisional Patent Application No. 62/983,421 entitled “Systems and Methods for Composing High Dimensional Fluorescence Microscopy Images In Subpixel Resolution” filed February 28, 2020, which is incorporated by reference herein in its entirety.

FIELD OF THE INVENTION

[0002] The present invention relates generally to image processing, including methods to compose high dimensional images with subpixel resolution. More particularly, the invention relates to computational methods for composing images from three-dimensional image data acquired by a microscope across multiple cycles of imaging.

BACKGROUND

[0003] There are a large number of problems with reproducibly generating microscopy images, especially at higher magnification and across time. When imaging in high magnification with small fields of view across time or repeated imaging cycles, microscope drift and variations in spatial positioning can prevent the proper assembly of a larger stitched image with high resolution. Assembling these images into biologically relevant images for downstream analysis is not trivial. This can often require the registration of images in 3D and algorithms to seamlessly stitch them together. Available algorithms for this are either slow, inaccurate, or degrade image quality.

SUMMARY OF THE INVENTION

[0004] Systems and methods for composing high dimensional microscopy images with subpixel resolution in accordance with various embodiments of the invention are disclosed. [0005] In one embodiment, a method for composing high dimensional microscopy images includes obtaining a plurality of microscopic images of a sample, where the images represent a three-dimensional stack of images representing a series of images that are tiled in each of X-, Y-, and Z-dimensions, determining a Z-position for each image in the plurality of microscopic images, determining best focus Z-positions for sub-regions within each image, determining a best focus surface based on the Z-positions of each image and the best focus Z-positions for each image, shifting at least some of the plurality of images into a central plane, deconvolving each stack of images, compensating forX/Y- drift of each tile based on the deconvolved images, registering at least some of the images in the plurality of images based on Z-positions and X/Y-drift for each image, and stitching at least some of the plurality of images to form an output image based on the registered images.

[0006] In a further embodiment, the plurality of microscopic images of a sample include a plurality of images representing light emission taken from a plurality of channels that represent band-filtered wavelengths of light, wherein the three-dimensional stack of images represents a series of images that are tiled in each of the X-, Y-, and Z- dimensions, wherein each tile position is imaged in each channel.

[0007] In another embodiment, the method further includes measuring misalignment in X-, Y-, and Z- planes between adjacent images in the plurality of images to determine X-, Y-, and Z-offsets, where the determining Z-positions step uses the measured misalignment to determine Z-positions.

[0008] In a still further embodiment, the measuring misalignment includes performing fast Fourier transformation on the plurality of images to obtain a frequency-domain representation of the images, determining a Hadamard product of each frequency-domain representation of the images, performing an inverse fast Fourier transformation to obtain a spatial correlation result, identifying a correlation peak between images in the plurality of images, and estimating a subpixel peak with three-dimensional (3D) quadratic interpolation.

[0009] In still another embodiment, determining Z-positions comprises optimizing Z- positions using iterative force convergence optimization. [00010] In a yet further embodiment, the iterative force convergence optimization is performed by initializing images positions to zero and summing forces represented by the X-, Y-, and Z-offsets.

[00011] In yet another embodiment, the plurality of microscopic images of a sample include a plurality of cycles, where each cycle represents a series of images that are tiled in each of the X-, Y-, and Z-dimensions.

[00012] In a further embodiment again, the method further includes measuring Z-drift across cycles.

[00013] In another embodiment again, measuring the Z-drift across cycles includes performing fast Fourier transformations on at least some of the plurality of images to obtain frequency-domain representations of the images, determining pairwise X-, Y-, and Z-offsets by measuring cross-correlation peaks between cycles, determining a Hadamard product of each of the frequency-domain representations of the images, perform inverse fast Fourier transformations to obtain a spatial correlation result, identifying a correlation peak between at least some of the images in the plurality of images, and estimating a subpixel peak with three-dimensional (3D) quadratic interpolation.

[00014] In a further additional embodiment, measuring the Z-drift across cycles further includes generating weighted X-, Y-, and Z-offsets for each cycle and shifting each cycle image stack into a position to correct for the measured Z-offset.

[00015] In another additional embodiment, determining a best focus surface is further based on the Z-drift across cycles.

[00016] In a still yet further embodiment, the plurality of microscopic images of a sample include a plurality of fluorescence images taken from a plurality of channels, where the three-dimensional stack of images represents a series of images that are tiled in each of the X-, Y-, and Z-dimensions, wherein each tile position is imaged in each channel during each cycle.

[00017] In still yet another embodiment, one channel in each cycle is used to image a specific structure in the sample.

[00018] In a still further embodiment again, determining a Z-position, determining a best focus surface, and compensating for X/Y-drift use one channel from each cycle that images the specific structure in the sample. [00019] In still another embodiment again, the method further includes applying Richardson-Lucy deconvolution to the images to reduce blurring.

[00020] In a still further additional embodiment, the method further includes reducing the number of Z-slices within the plurality of images to remove Z-slices that do not contain the sample.

[00021] In still another additional embodiment, the method further includes performing an extended depth of field calculation on at least some of the plurality of images.

[00022] In a yet further embodiment again, the method further includes compressing high dynamic range image data.

[00023] In yet another embodiment again, a non-transitory machine-readable media containing processor instructions, where execution of the instructions causes a processor to perform a process to compose high dimensional microscopy images including obtaining a plurality of microscopic images of a sample, where the images represent a three-dimensional stack of images representing a series of images that are tiled in each of X-, Y-, and Z-dimensions, determining a Z-position for each image in the plurality of microscopic images, determining best focus Z-positions for sub-regions within each image, determining a best focus surface based on the Z-positions of each image and the best focus Z-positions for each image, shifting at least some of the plurality of images into a central plane, deconvolving each stack of images, compensating for X/Y-drift of each tile based on the deconvolved images, registering at least some of the images in the plurality of images based on Z-positions and X/Y-drift for each image, and stitching at least some of the plurality of images to form an output image based on the registered images.

[00024] In a yet further additional embodiment, a system for composing high dimensional microscopy images includes a processor and a memory readable by the processor, where the memory contains instructions that when read by the processor direct the processor to obtain a plurality of microscopic images of a sample, wherein the images represent a three-dimensional stack of images representing a series of images that are tiled in each of X-, Y-, and Z-dimensions, determine a Z-position for each image in the plurality of microscopic images, determine best focus Z-positions for each image, determine a best focus surface based on the Z-positions of each tile and the best focus Z-positions for each tile, shift at least some of the plurality of images into a central plane, deconvolve each image, compensate for X/Y-drift of each tile based on the deconvolved images, register at least some of the plurality of images based on Z-positions and X/Y- drift for each image, and stitch at least some of the plurality of images to form an output image based on the registered images.

[00025] Additional embodiments and features are set forth in part in the description that follows, and in part will become apparent to those skilled in the art upon examination of the specification or may be learned by the practice of the disclosure. A further understanding of the nature and advantages of the present disclosure may be realized by reference to the remaining portions of the specification and the drawings, which forms a part of this disclosure.

BRIEF DESCRIPTION OF THE DRAWINGS

[00026] These and other features and advantages of the present invention will be better understood by reference to the following detailed description when considered in conjunction with the accompanying drawings where:

[00027] Figure 1 illustrates a method for composing high dimensional images in accordance with various embodiments.

[00028] Figure 2 illustrates a multi-cycle imaging process in accordance with various embodiments.

[00029] Figure 3 illustrates three-dimensional offsets between image tiles in accordance with various embodiments.

[00030] Figure 4 illustrates an exemplary force convergence model in accordance with various embodiments.

[00031] Figure 5 illustrates inter-cycle Z-offsets between image tiles in accordance with various embodiments.

[00032] Figure 6 illustrates surface fitting to imaged tissue in accordance with various embodiments.

[00033] Figures 7A-7B illustrate image deconvolution in accordance with various embodiments. [00034] Figures 8A-8C illustrate chromatic aberration correction in accordance with various embodiments.

[00035] Figures 9A-9B illustrate image stitching in accordance with various embodiments.

[00036] Figures 10A-10D illustrate exemplary image reconstruction in accordance with various embodiments.

DETAILED DESCRIPTION

[00037] The embodiments of the invention described herein are not intended to be exhaustive or to limit the invention to precise forms disclosed. Rather, the embodiments selected for description have been chosen to enable one skilled in the art to practice the invention.

[00038] Microscopy is a valuable tool in studying tissues and cells of an organism. Flowever, microscopes are able to focus on a single plane at a time, when tissues can be three dimensional, thus disallowing imaging of certain features within a single view. Even when imaging across multiple planes, diffuse light from out of focus and/or background elements can reduce resolution in any resulting image. Many embodiments are capable of correcting images across time or multiple imaging cycles, such as when imaging tissue with multiple probes across multiple cycles of imaging (e.g., CODEX, CO-Detection by indEXing). ( See e.g., Goltsev, Y. et al. Deep Profiling of Mouse Splenic Architecture with CODEX Multiplexed Imaging. Cell 174, 968-981. e15 (2018); the disclosure of which is hereby incorporated by reference in its entirety.) During multi-cycle imaging, many problems become more complex as imaged tissue can move slightly between cycles, making stitching even more difficult than during a single imaging cycle.

[00039] Turning now to the drawings, many embodiments are directed to methods to apply corrections to images to allow for efficient and highly accurate registration and stitching. Numerous embodiments are capable of applying image corrections to any multichannel imaging system possessing a reference channel, including single-cycle and multi-cycle imaging. Various embodiments can be applied to Raman spectrometry, light and brightfield microscopy, spatial mass spectrometry imaging, multispectral satellite imagery, and other relevant forms of imaging known in the art. In certain embodiments, the images are obtained from band-filtered wavelengths of light, such as fluorescence microscopy. A number of embodiments are capable of providing sub-pixel resolution while shifting images during processing, which provides higher-precision alignment across image dimensions, better maintenance of image integrity and quality, and sub-frame shifting of images to adjust images to a single in-focus plane. In several embodiments, alignment is performed by transforming images into the Fourier or frequency domain and determining the offset(s) that yields the greatest correlation between transformed images. Additional embodiments are able to have faster implementation by using graphics processing units (GPUs) over central processing units (CPUs).

[00040] Turning to Figure 1 , a method 100 for composing high dimensional images with subpixel resolution in accordance with an embodiment of the invention is illustrated. Images are obtained (102). A number of embodiments obtain multiple images tiled from at least one dimension (e.g., images tiled in the X-dimension). Additional embodiments obtain multiple images tiled from at least two dimension (e.g., images tiled in both X- and Y-dimensions). Further embodiments obtain multiple images that are tiled from at least 3 dimensions (e.g., images tiled across X-, Y-, and Z-dimensions). In some embodiments, the images that are obtained have overlap between the images, such that there could be 10%, 20%, 30%, 40%, 50%, or more overlap between images that are obtained by moving the microscope in any dimension (e.g., adjacent images in the X-dimension has 30% overlap with each other).

[00041] In some embodiments, the images are obtained directly from a microscope. In several embodiments, images are obtained from a storage device, such as a hard drive or server. Microscopic images in numerous embodiments are obtained from wide field microscopy, while certain embodiments obtain microscopic images using confocal microscopy. In further embodiments, the microscopic images are obtained using fluorescence microscopy. A number of embodiments can obtain microscopic images from a series of cycles of imaging, where one cycle images a specimen (e.g., tissue) across one or more channels, where a channel represents a specified range of light wavelengths (e.g., red channel (~620-750nm), green channel (~495-570nm), blue channel (~450- 495nm), etc.). Certain embodiments can perform imaging across at least 5, 10, 15, 20, or more cycles, with each cycle possessing imaging data from at least 1 , 2, 3, 4, or more channels. An example of multi-cycle imaging is illustrated in Figure 2, which shows a single field of view imaged across 6 cycles. Further, a constant reference channel is used across each cycle as DAPI used to stain nuclei in each cell, while other channels vary with different probes at each cycle.

[00042] Returning to Figure 1 , in many embodiments, the misalignment between adjacent images (e.g., tiles) in a z-stack is measured (104). In a number of embodiments image data from a single channel (e.g., blue channel, green channel, red channel, etc.) in a single cycle is used to measure misalignment in the obtained images. The images from the single channel in a single cycle can be dark and/or flatfield corrected. In some embodiments, overlapping edges between adjacent images can be cross-correlated to measure a peak to yield an offset (e.g., DC, DU, DZ offsets). Figure 3 illustrates examples of X-, Y-, and Z- offsets from neighboring images within a single cycle. In several embodiments, offsets further exist for a particular position (e.g., XYZ coordinate) between cycles (e.g., offsets between cycle 2 and cycle 3).

[00043] In a number of embodiments, the cross-correlation is measured by performing Fourier transformation (such as fast Fourier transformation) to obtain a frequency-domain representation of the images. In a number of embodiments, a Fladamard product off^* and g are determined for each frequency-domain image, where f^* represents the complex conjugate of f. In certain embodiments, a filtering function is applied to discard or attenuate high- or low- frequency components in the frequency-domain representations. As can readily be appreciated, any of a variety of filters can be applied as appropriate to the requirements of specific applications. Certain embodiments can also perform an inverse Fourier transformation (including inverse fast Fourier transformation) to obtain a spatial correlation result. Once the spatial correlation result is determined, various embodiments can identify a correlation peak and an associated subpixel shift using three- dimensional (3D) quadratic interpolation. The X-, Y-, and Z- shifts and correlation scores for each tile relative to its neighbors can be stored as a measure of drift. In several embodiments the correlation scores are multiplied by a function of the image means and standard deviations. [00044] Turning back to Figure 1 , in many embodiments, tile Z-positions are determined (106) based on a consensus of Z-offsets in a single cycle in order to reconcile tile positions with confidence from previously determined offset measurements (104). In certain embodiments, measurements that are outliers for the X-, Y-, and Z-offsets and/or correlation scores are discarded. In a number of embodiments, the tile positions are optimized via iterative force convergence. In a number of embodiments, the iterative force convergence optimization is performed by initializing tile positions to zero and summing forces represented by the X-, Y-, and Z-offsets. Each tile position can be updated proportionally to a weighted average force for each tile in a number of embodiments, and scaling factors are decreased per iteration to achieve convergence of tiles into a reconciled position. The updated tile positions can then be recorded or saved. As can readily be appreciated, any of a variety of processes can be utilized to determine the relative positions of tiles based upon offset correlation measurements as appropriate to the requirements of specific applications in accordance with various embodiments of the invention. Figure 4 illustrates an exemplary demonstration of an iterative force convergence model to determine and/or optimize Z-positions based on Z-offsets.

[00045] Referring back to Figure 1 , in several embodiments, z-drift is measured (108) across cycles. For example, Figure 5 illustrates a Z-offset that occurs between imaging cycles (e.g., cycle m and cycle n). In certain embodiments, the images from a single channel across multiple cycles can be dark and/or flatfield corrected. In a number of embodiments, the across cycle drift is measured by performing Fourier transformation (such as fast Fourier transformation) to obtain a frequency-domain representation of the images. In some embodiments, image stacks (e.g., z-stacked images) are shifted to the correct identified (106) z-positions. The frequency-domain representation can be multiplied by a shifting kernel computed form fractional pixel Z-offsets. In some embodiments z-stacked images are shifted in the frequency domain by kernel(z) = e^A(k ^* 2pi ^* (z ^* dz / number_of_slices). Once image stacks are corrected, a number of embodiments perform comparisons between cycle pairs to determine pairwise X-, Y-, and Z-offsets by measuring cross-correlation peaks between cycles. In a number of embodiments, a Fladamard product of f^* and g are determined for each image, where f^* represents the complex conjugate of f. In certain embodiments, the filtering function is introduced to discard or attenuate high- or low- frequency components in the frequency domain representations. As can readily be appreciated, any of a variety of filters can be applied as appropriate to the requirements of specific applications. In certain embodiments, inverse Fourier transformation (including inverse fast Fourier transformation) are performed to obtain a spatial correlation result. Once the spatial correlation result is determined, various embodiments can find a correlation peak and estimate a subpixel peak with three-dimensional (3D) quadratic interpolation. The X-, Y-, and Z- shifts and correlation scores for each tile relative to its neighbors can be stored as a measure of drift. In further embodiments, the correlation scores can be multiplied by a function of the image means and standard deviations. Several embodiments can also determine a “best” cycle, which has the greatest total correlation weight. Additional embodiments compare pairwise offsets with the best cycle and reject any offsets that disagree with the best cycle. Finally, certain embodiments generate weighted X-, Y-, and Z-offsets for each cycle. Once weighted offsets are determined, some embodiments will shift each cycle image stack into a position to correct for the measured Z-offset between cycles. Further embodiments also perform an inverse Fourier transformation (including inverse fast Fourier transformation) to obtain spatial-domain versions of the position- and cycle-corrected image stacks.

[00046] Several embodiments determine (110) a best focus slice (Z-position) for each tile stack. In many embodiments, image stacks are diced into sub-tile stacks representing a subset of pixels of each frame. Some embodiments dice the image stacks into 256X256 pixels in the X- and Y-dimensions. As can readily be appreciated the specific dimensions and corresponding number of diced image stacks are largely determined by the requirements of a particular application. In certain embodiments, each sub-tile position is compared at each Z-position in a pairwise manner across imaging cycles to compute a focus score for a particular Z-position by cross-correlation. To determine best focus, certain embodiments perform Fourier transformation (such as fast Fourier transformation) to obtain a frequency-domain representation of the images. In some embodiments, pairs of cycles are compared to determine relative X-, Y-, and Z-offsets by measuring cross correlation peaks between cycles. In a number of embodiments, a Fladamard product of f^* and g are determined for each image, where f^* represents the complex conjugate of f. Certain embodiments will further perform inverse Fourier transformation (including inverse fast Fourier transformation) to obtain a spatial correlation result. Once determining the spatial correlation result, various embodiments will find a correlation peak and estimate a subpixel peak with two-dimensional (2D) quadratic interpolation. The X-, Y-, and Z- shifts and correlation scores for each tile relative to its neighbors can be stored as a measure of drift in many embodiments. In further embodiments the correlation scores are multiplied by a function of the image means and standard deviations. Additional embodiments can accumulate correlation scores for each sub-tile position and Z-position across all pairs of cycles. The total score of for each sub-tile position and each Z-position can be recorded or saved in many embodiments. Although specific processes are described above for determining a best focus, any of a variety of processes can be utilized to determine best focus as appropriate to the requirements of specific applications in accordance with various embodiments of the invention.

[00047] In certain embodiments, a best focus surface is determined (112) in conjunction with corresponding Z-offsets for an imaged region. Many embodiments utilize Z-offsets, cross-correlation information, and best focus scores, such as those generated using processes similar to those outlined above, and a peak fitting algorithm, such as (but not limited to) quadratic interpolation, to determine the best focus slice (e.g., Z-position) for each sub-tile position. In a number of embodiments, the best focus surface is identified by fitting a smooth 2D surface to the positional best focus slice data. In certain embodiments, sub-tile fitted best focus surface offsets are combined with cycle drift (e.g., Z-drift) and the resulting offsets are saved or recorded. An exemplary focused surface is illustrated in Figure 6, where surface 602 has been fitted to tissue 604.

[00048] Returning to Figure 1 , many embodiments deconvolve and Z-shift images (114) into a central plane. In many embodiments, the process looks at an image stack (e.g., Z- stack) from a single position, cycle, and channel at a time. Images can be dark and/or flatfield corrected. In many embodiments, image stacks are diced into sub-tile stacks representing a subset of pixels. Some embodiments will dice the image stacks into 256X256 pixels in the X- and Y-dimensions. As can readily be appreciated, the specific dicing process is largely determined based upon the requirements of a given application. In a number of embodiments, a Fourier transformation (such as fast Fourier transformation) is performed on the images to obtain a frequency-domain representation of the images. In some embodiments, each sub-tile image stack is shifted to move the best focus slice into the middle of the sub-tile stack using the previously determined positional Z-shift (112). The frequency-domain representation can be multiplied by a shifting kernel computed form fractional pixel Z-offsets. In some embodiments z-stacked images are shifted in the frequency domain by kernel(z) = e^A(k ^* 2pi ^* (z ^* dz / number_of_slices). Further embodiments also apply a de-ringing kernel of a specified strength to reduce Gibbs artifacts which are a natural consequence applying frequency- domain shifting with fractional pixel shifts. As can readily be appreciated, any of a variety of image processing techniques can be utilized to remove artifacts and/or improve image quality as appropriate to the requirements of specific applications.

[00049] In a number of embodiments, Richardson-Lucy deconvolution is applied (116) to the corrected images. Deconvolution is applied in such embodiments to reduce blurring within images. Additional embodiments generate point spread functions for the deconvolution using a Gibson-Lanni, Born-Wolf, or another model. Additional embodiments perform blind deconvolution by updating the point spread function with each iteration of deconvolution. Further embodiments utilize a “semi-blind” deconvolution, where the point spread function starts with a “best guess” prior to updating with iterations of deconvolution. Figures 7A-7B illustrate exemplary effects of deconvolution in accordance with various embodiments. Specifically, Figure 7AS illustrates an initial image prior to deconvolution, and Figure 7B illustrates the same image, once deconvolution has been applied in accordance with various embodiments. Figure 7B thus shows that many embodiments are capable of reducing blur in collected images.

[00050] In many embodiments, certain Z-slices will not contain tissue or the sample of interest, thus these “empty” Z-slices can be removed, as they do not contain relevant information. Thus, additional embodiments trim deconvolved image stacks to reduce the number of output Z-slices(118). Further embodiments perform an extended depth of field (EDF) calculation (120). The specific process utilized to achieve EDF is largely dependent upon the requirements of a given application. ( See e.g., Forster, et al. “Complex wavelets for extended depth-of-field: A new method for the fusion of multichannel microscopy images,” Microsc. Res. Tech. 65:33-42, 2004; the disclosure of which is incorporated by reference herein in its entirety.) Deconvolution in accordance with various embodiments, allows tissue and fitted surface to be shifted in three-dimensional space, such that the tissue that may be on an incline within the slide (e.g., Figure 6) can be stored as a more planar or flattened image relative to the slide. In addition to making the resultant image construction more convenient as having a surface to be more planar, such shifting can also reduce storage needs for the resulting image construct by reducing the range of data necessary to store the deconvolved image.

[00051] In several embodiments, floating-point image data can be compressed (122) into a 16-bit format. To preserve the informational content from 32-bit deconvolved images, certain embodiments identify intensity values as constants such as “low,” “high,” and “A,” and negative values can be stored as 0 values. In various embodiments, “low” is set as an inflection point where values lower than “low” are stored linearly and values higher than “high” are stored as 65534, whereas values between “low” and “high” are transformed by a predetermined function that compresses the range of image values between “low” and “high” into the output range of “low” to 65533. Certain embodiments will use the function (LOW + 1 + A ^* (log(V - LOW - 0.5 - 1/(exp(1/A)-1 )) - log(1/(exp(1/A)- 1 ))) for a chosen value of “A” to store intensity values between “low” and “high”. Further embodiments will store NaN values as 65535. As can readily be appreciated, any of a variety of compression functions can be utilized as appropriate to the requirements of a given application in accordance with various embodiments of the invention.

[00052] Several embodiments compensate for X/Y-drift (124). Images from a single channel and position across cycles can be utilized to compensate for X/Y drift. A number of embodiments perform a Fourier transformation (such as fast Fourier transformation) to obtain a frequency-domain representation of the images. For each integer slice offset in a range (e.g., maximum slice offset in the negative and positive direction), certain embodiments compare pairs of cycles at each Z-slice (with slice offset) to determine a pairwise X/Y-offset by measuring a cross-correlation peak. The cross-correlation peak can be determined computing a Hadamard product of frequency-domain images f^* and g, where f^* represents the complex conjugate of f. In certain embodiments, a filtering function is introduced to discard or attenuate high- or low- frequency components in the frequency-domain representations. Certain embodiments will further perform an inverse Fourier transformation (including inverse fast Fourier transformation) to obtain a spatial correlation result. Once the spatial correlation result is determined, various embodiments find a correlation peak and estimate a subpixel peak with 2D quadratic interpolation. The X-, and Y- shifts and correlation scores for each cycle pair with are computed for each slice offset. Many embodiments determine the best slice offset for each cycle pair as the greatest total correlation score. More embodiments determine a “best” cycle with the greatest total correlation score. The best cycle and each pairwise offset can be compared to reject any offsets that disagree with the best offset. Finally, the X-, Y-, and Z-offsets for each cycle can be determined and recorded or stored.

[00053] Image registration can then be performed (126). In many embodiments, drift compensation information is used to determine X-, Y-, and Z-drifts between cycles, such as previously determined drift information. Any outlier (e.g., low weight or large offsets) can be discarded. Images corresponding to a particular Z-slice can be selected and for each grid position, and overlapping portions of images (e.g., overlapping in the X- and/or Y-dimension) can be loaded for each channel. Further embodiments can register the loaded images by cross-correlation. During cross-correlation, many embodiments perform a Fourier transformation (such as fast Fourier transformation) to obtain a frequency-domain representation of the images. A Fladamard product of f^* and g can be determined for each frequency-domain image pair, where f^* represents the complex conjugate of f. In certain embodiments, a filtering function is applied to discard or attenuate high- or low- frequency components in the frequency domain representations. As can readily be appreciated, any of a variety of filtering processes can be performed as appropriate to the requirements of a given application. In several embodiments, an inverse Fourier transformation (including inverse fast Fourier transformation) is performed to obtain a spatial correlation result. Further embodiments can multiply the spatial correlation matrix by a gaussian 'bias matrix' centered around the expected position. In a number of embodiments, a correlation peak is identified by estimating a subpixel peak with 2D quadratic interpolation. The correlation peak can also be refined using processes including (but not limited to) discrete Fourier transform upsampling. ( See e.g., Guizar- Sicairos, et al. “Efficient subpixel image registration algorithms,” Opt Lett. 2008 Jan 15;33(2): 156-8; the disclosure of which is incorporated herein by reference in its entirety.) [00054] Further embodiments combine X- and Y-offsets for each cycle and channel after adding chromatic aberration terms for each channel. For example, Figure 8A illustrates how different wavelengths of light refract differently through lenses or other optics, causing differences in the camera sensor F. As such, Figure 8B illustrates how imaging of the same point in three different channels can create aberration based on the refractive differences. Certain embodiments determine residual chromatic aberration by accumulating channel-channel offset differences across all image registrations to create a corrected image, such as illustrated in the exemplary image of Figure 8C.

[00055] Additional embodiments can store or record weighted X- and Y-shifts and correlation scores for overlapping edges. Certain embodiments determine global X/Y shifts for the X/Y-plane as weighted averages of all measurements. Tile positions can be set based on this global alignment. In a number of embodiments, the tile positions are optimized via an iterative force convergence algorithm. Further embodiments can iteratively refine tile alignment to obtain better image registration. Once image registration is complete, tile positions can be stored or recorded. Further embodiments measure a registration error and produce an error score.

[00056] Referring back to Figure 1 , many embodiments perform high dimensional image stitching (128). For example, as illustrated in Figures 9A-9B, image stitching are illustrated. Specifically, Figure 9A illustrates an erroneous stitching effect (pages.nist.gov/MIST) having an error that creates a vertical line between neighboring tiles (see arrows). Many embodiments are capable of stitching the images with much better precision that does not have similar artifacts due to subpixel precision, such as illustrated in the example of Figure 9B. The tile registration information can be utilized for each channel and each cycle, as appropriate. For each tile, image data for a correct Z- slice at that position can be loaded, in many embodiments. Chromatic aberration correction can be applied to the image. The image can be shifted to its fractional pixel position by performing a Fourier transformation (such as fast Fourier transformation) to obtain a frequency-domain representation of the images. The frequency-domain representation can be multiplied by a shifting kernel computed form fractional pixel X- and Y-offsets. Additional embodiments can perform frequency domain shifting by kernel(x,y) = e^A(i ^* 2pi ^* (y ^* dy/h + x ^* dx/w)). In certain embodiments, resulting shifted representations are multiplied by a filtering function to attenuate high-frequency components. Further embodiments can apply a de-ringing function of a specified strength to reduce Gibbs artifacts, which are a consequence of applying frequency-domain shifting with fractional shifts. Additional embodiments apply an artifact-rejecting compositing kernel to produce a crisp and clean output image by selectively blending filtered and unfiltered versions of the shifted image. As can readily be appreciated, any of a variety of post-processing processes can be applied to remove artifacts and/or increase image quality as appropriate to the requirements of specification applications in accordance with various embodiments of the invention. Once the above image shifting is completed, many embodiments place the shifted image in a virtual canvas at the correct position. Additional embodiments blend overlapping areas by interpolating overlapping regions using a function such as (but not limited to) a ramp function. Additional embodiments can perform background subtraction of the virtual canvas. In some embodiments, tile registration information can be calculated for one ‘reference’ z-slice and applied to other z-slices for image stitching.

[00057] Once images are stitched, various embodiments utilize one or more machine learning models to annotate information from the sample based on the information regarding each imaged channel. For example, if each channel imaged represents a unique probe, biomarker, antibody, etc., such information can identify individual cells based on localization of such probes and combinations of probes on each cell. These annotations can be used for further applications of high dimensional imaging in accordance with various embodiments.

[00058] It should be noted that a variety of embodiments perform various operations multiple times, where appropriate. Where steps are relevant to a single cycle, these steps can be reiterated for each cycle. Additionally, certain steps may be omitted. For example, if applying certain embodiments to images coming from a single cycle, determining offsets between cycles would not be determinable. Further, certain embodiments can perform various processes described above with reference to Figure 1 steps in a different order and/or certain steps are omitted and/or performed simultaneously or in an overlapping fashion. [00059] In certain embodiments, a non-transitory machine-readable media is provided containing processor instructions, where execution of the instructions causes a processor to execute a method similar to any of the methods described above with reference to Figure 1. Further, image processing systems in accordance with many embodiments of the invention include a processor and memory, where the memory contains software that configures the processor to perform a process similar to one of the processes described above with reference to Figure 1 .

[00060] Many embodiments are capable of providing subpixel precision in aligning and constructing of images after multiple cycles of microscopic imaging. Figures 10A-10D illustrate exemplary alignments of images obtained from multi-cycle imaging. Specifically, Figures 10A-10B illustrate alignments of images from a second and fifteenth cycle of multi-cycle imaging, while Figures 10C-10D illustrate exemplary outputs of certain embodiments of image reconstruction. Specifically, Figure 10C illustrates exemplary output of an embodiment of the cycles aligned in Figure 10A, and Figure 10D illustrates exemplary output of an embodiment of the cycles aligned in Figure 10B.

DOCTRINE OF EQUIVALENTS

[00061] Flaving described several embodiments, it will be recognized by those skilled in the art that various modifications, alternative constructions, and equivalents may be used without departing from the spirit of the invention. Additionally, a number of well-known processes and elements have not been described in order to avoid unnecessarily obscuring the present invention. Accordingly, the above description should not be taken as limiting the scope of the invention.

[00062] Those skilled in the art will appreciate that the presently disclosed embodiments teach by way of example and not by limitation. Therefore, the matter contained in the above description or shown in the accompanying drawings should be interpreted as illustrative and not in a limiting sense. The following claims are intended to cover all generic and specific features described herein, as well as all statements of the scope of the present method and system, which, as a matter of language, might be said to fall therebetween.

Claims

What is claimed is:

1. A method for composing high dimensional microscopy images comprising: obtaining a plurality of microscopic images of a sample, wherein the images represent a three-dimensional stack of images representing a series of images that are tiled in each of X-, Y-, and Z-dimensions; determining a Z-position for each image in the plurality of microscopic images; determining best focus Z-positions for sub-regions within each image; determining a best focus surface based on the Z-positions of each image and the best focus Z-positions for each image; shifting at least some of the plurality of images into a central plane; deconvolving each stack of images; compensating forX/Y-drift of each tile based on the deconvolved images; registering at least some of the images in the plurality of images based on Z- positions and X/Y-drift for each image; and stitching at least some of the plurality of images to form an output image based on the registered images.

2. The method of claim 1, wherein the plurality of microscopic images of a sample include a plurality of images representing light emission taken from a plurality of channels that represent band-filtered wavelengths of light, wherein the three- dimensional stack of images represents a series of images that are tiled in each of the X-, Y-, and Z-dimensions, wherein each tile position is imaged in each channel.

3. The method of claim 1, further comprising measuring misalignment in X-, Y-, and Z- planes between adjacent images in the plurality of images to determine X-, Y-, and Z-offsets, wherein the determining Z-positions step uses the measured misalignment to determine Z-positions.

4. The method of claim 3, wherein the measuring misalignment step comprises the steps of: performing fast Fourier transformation on the plurality of images to obtain a frequency-domain representation of the images; determining a Hadamard product of each frequency-domain representation of the images; performing an inverse fast Fourier transformation to obtain a spatial correlation result; identifying a correlation peak between images in the plurality of images; and estimating a subpixel peak with three-dimensional (3D) quadratic interpolation.

5. The method of claim 1, wherein determining Z-positions comprises optimizing Z- positions using iterative force convergence optimization.

6. The method of claim 5, wherein the iterative force convergence optimization is performed by initializing images positions to zero and summing forces represented by the X-, Y-, and Z-offsets.

7. The method of claim 1, wherein the plurality of microscopic images of a sample include a plurality of cycles, wherein each cycle represents a series of images that are tiled in each of the X-, Y-, and Z-dimensions.

8. The method of claim 7, further comprising measuring Z-drift across cycles.

9. The method of claim 8, wherein the measuring the Z-drift across cycles comprises: performing fast Fourier transformations on at least some of the plurality of images to obtain frequency-domain representations of the images; determining pairwise X-, Y-, and Z-offsets by measuring cross-correlation peaks between cycles; determining a Fladamard product of each of the frequency-domain representations of the images; perform inverse fast Fourier transformations to obtain a spatial correlation result; identifying a correlation peak between at least some of the images in the plurality of images; and estimating a subpixel peak with three-dimensional (3D) quadratic interpolation.

10. The method of claim 9, wherein the measuring the Z-drift across cycles further comprises: generating weighted X-, Y-, and Z-offsets for each cycle; and shifting each cycle image stack into a position to correct for the measured Z- offset.

11. The method of claim 8, wherein the determining a best focus surface is further based on the Z-drift across cycles.

12. The method of claim 7, wherein the plurality of microscopic images of a sample include a plurality of fluorescence images taken from a plurality of channels, wherein the three-dimensional stack of images represents a series of images that are tiled in each of the X-, Y-, and Z-dimensions, wherein each tile position is imaged in each channel during each cycle.

13. The method of claim 12, wherein one channel in each cycle is used to image a specific structure in the sample.

14. The method of claim 13, wherein determining a Z-position, determining a best focus surface, and compensating for X/Y-drift use one channel from each cycle that images the specific structure in the sample.

15. The method of claim 1, further comprising applying Richardson-Lucy deconvolution to the images to reduce blurring.

16. The method of claim 15, further comprising reducing the number of Z-slices within the plurality of images to remove Z-slices that do not contain the sample.

17. The method of claim 1, further comprising performing an extended depth of field calculation on at least some of the plurality of images.

18. The method of claim 1 , further comprising compressing high dynamic range image data.

19. A non-transitory machine-readable media containing processor instructions, where execution of the instructions causes a processor to perform a process to compose high dimensional microscopy images comprising: obtaining a plurality of microscopic images of a sample, wherein the images represent a three-dimensional stack of images representing a series of images that are tiled in each of X-, Y-, and Z-dimensions; determining a Z-position for each image in the plurality of microscopic images; determining best focus Z-positions for each image; determining a best focus surface based on the Z-positions of each tile and the best focus Z-positions for each tile; shifting at least some of the plurality of images into a central plane; deconvolving each image stack; compensating forX/Y-drift of each tile based on the deconvolved images; registering at least some of the plurality of images based on Z-positions and X/Y-drift for each image; and stitching at least some of the plurality of images to form an output image based on the registered images.

20. A system for composing high dimensional microscopy images comprising: a processor; and a memory readable by the processor, wherein the memory contains instructions that when read by the processor direct the processor to: obtain a plurality of microscopic images of a sample, wherein the images represent a three-dimensional stack of images representing a series of images that are tiled in each of X-, Y-, and Z-dimensions; determine a Z-position for each image in the plurality of microscopic images; determine best focus Z-positions for each image; determine a best focus surface based on the Z-positions of each tile and the best focus Z-positions for each tile; shift at least some of the plurality of images into a central plane; deconvolve each image; compensate forX/Y-drift of each tile based on the deconvolved images; register at least some of the plurality of images based on Z-positions and X/Y-drift for each image; and stitch at least some of the plurality of images to form an output image based on the registered images.