WO2010105015A2

WO2010105015A2 - Methods, systems, and computer readable media for microscopy tracking

Info

Publication number: WO2010105015A2
Application number: PCT/US2010/026913
Authority: WO
Inventors: Brian Eastwood; Russell M. Ii Taylor; Richard Superfine; Lamar Mair; Marc Niethammer
Original assignee: The University Of North Carolina At Chapel Hill
Priority date: 2009-03-11
Filing date: 2010-03-11
Publication date: 2010-09-16
Also published as: WO2010105015A3

Abstract

Methods and systems for microscopy tracking are disclosed. According to one aspect of the subject matter described herein, a system for microscopy tracking is disclosed. The system can comprise a specimen holder or covering member for holding or covering a specimen being imaged using a microscope and a pattern formed in or on the specimen holder or covering member. The pattern may be structured for tracking locations of features of the specimen in images of the specimen obtained from the microscope.

Description

DESCRIPTION

METHODS, SYSTEMS, AND COMPUTER READABLE MEDIA FOR

MICROSCOPY TRACKING

PRIORITY CLAIM

This application claims the benefit of U.S. Provisional Patent Application Serial No. 61/159,330, filed March 11 , 2009, the disclosure of which is incorporated herein by reference in its entirety.

GOVERNMENT INTEREST

This work was supported at least in part by a grant from the National Institutes of Health (Grant No. P41-EB002025-23 to 25A1). Thus, the U.S. Government may have certain rights in the presently disclosed subject matter.

TECHNICAL FIELD

The subject matter described herein relates to microscopy. More specifically, the subject matter relates to methods, systems, and computer readable media for microscopy tracking.

BACKGROUND

In general, traditional microscopes require researchers to forfeit high magnification for field of view (FOV). For example, using a microscope at high magnification is excellent at enabling a user to view fine structures in small areas of a specimen, however, viewing large areas at high magnification requires the user to continuously move a microscope stage. Microscopists are often interested in both of these aspects simultaneously, for example to understand muscle morphology. A single view of large areas of the specimen at high magnification may be built by fusing multiple high magnification images of small areas together to form an image mosaic. This is analogous to creating a panoramic view by stitching multiple small images together in macroscopic photography. Stitching multiple images together requires registering each constituent image to a single reference frame. One may obtain this information by either tracking the imaging device (i.e., the stage in microscopy) or tracking image based features.

Current methods for image based tracking include using electronic/motor driven stages or using contrast within an image to track and recognize feature points. Using contrast within an image often requires assuming a specimen is fixed. However, each method presents problems. For example, not all microscopes are equipped with electronic/motor driven stages, as such stages can be very expensive. Even if purchased, high quality stages may be inaccurate, misaligned, and susceptible to drift. Tracking methods depending on sufficient contrast within an image to track or recognize feature points can also be challenging as live specimens may move or specimens may contain sparse areas of low contrast. Thus extracting stage motion from specimen motion and tracking features in large areas of low contrast can be challenging. For example, if a specimen moves, such as tracking mucus transport across cell cultures or investigating cell motility, the specimen can move further than one field of view will allow observation. As a result, it becomes more difficult to separate specimen motion from stage motion.

For at least these reasons, a need for methods and systems for microscopy tracking exists.

SUMMARY

According to another aspect of the subject matter described herein, a method for microscopy tracking is disclosed. The method can comprise providing a specimen holder or covering member for holding or covering a specimen being imaged using a microscope, tracking, using the pattern, locations of features of the specimen in images of the specimen obtained from the microscope. The pattern may be located in or on the sample holder or covering member, and may be structured to allow tracking of locations of features of the specimen.

The subject matter described herein for microscopy tracking may be implemented using a non-transitory computer readable medium having stored thereon executable instructions that when executed by the processor of a computer control the computer to perform steps. Exemplary computer readable media suitable for implementing the subject matter described herein include chip memory devices, disk memory devices, programmable logic devices, and application specific integrated circuits. In addition, a computer readable medium that implements the subject matter described herein may be located or implemented on a single device or computing platform or may distributed across plural devices or computing platforms.

BRIEF DESCRIPTION OF THE DRAWINGS

The subject matter described herein will now be explained with reference to the accompanying drawings of which:

Figure 1 is a block diagram of an exemplary system for microscopy tracking according to an embodiment of the subject matter described herein; Figure 2 is a flow chart of an exemplary process for microscopy tracking according to an embodiment of the subject matter described herein; Figure 3A is a micrograph of a square grid pattern; Figure 3B is the Fourier transform of the square grid pattern shown in Figure 3A; Figure 4A is micrograph of a square grid pattern;

Figure 4B is an illustration of a subset of Fourier transform magnitude for the pattern in Figure 4A;

Figure 4C is a graph of the magnitude along axis U; Figure 4D is a graph of the magnitude along axis V; Figure 5A is a graph of the Discrete Fourier transform of a single-frequency sinusoid;

Figure 5B is a graph of windowing effect on the discrete signal from Figure 5A; Figure 5C is a graph of the Discrete Fourier transform of a single- component signal;

Figure 5D is a graph of windowing effect on the discrete signal from Figure 5C; Figure 6A is a graph of an original input signal;

Figure 6B is a graph of the partial sum of the signal in Figure 6A;

Figure 6C is a graph showing the number of contributions for each element of Figure 6B;

Figure 6D is a graph of the partial sum of Figure 6B reweighted by Figure 6C;

Figure 6E is a graph of the Fourier basis function with target wavelength;

Figure 6F is a graph of the weighted signal 6D multiplied by Fourier basis function 6E;

Figure 7A is a micrograph of the axis-aligned input image with target wavelength marked by dashed vertical lines;

Figure 7B is a micrograph showing a partial sum of the micrograph in Figure 7A by target wavelength;

Figure 7C is a micrograph showing the number of contributions for each element of 7B; Figure 7D is a micrograph showing the partial sum of Figure 7B reweighted by Figure 7C;

Figure 7E is a micrograph of the Fourier basis function with target wavelength;

Figure 7F is an illustration of the weighted signal of Figure 7D multiplied by the Fourier basis function of Figure 7E;

Figure 8A is a graph of the mean phase error for all methods using a rectangular window discrete Fourier transform (DFT);

Figure 8B is a graph of the mean phase error from Figure 8A for the weighted phase correlation (WPC) method only; Figure 9A is a graph of the mean phase error for a corrupted signal using a rectangular window DFT;

Figure 9B is a graph of the mean phase error from Figure 9A, for the WPC method only; Figure 10A is an image of a grid predicted by geometrical optics;

Figure 10B is an image of a grid predicted by Fourier optics;

Figure 10C is an image simulated using quantized to whole-integer intensity values; Figure 11A is a graph of mean displacement error for different noise models using WPC;

Figure 11B is a graph of mean displacement error for different noise models using model-based spatial correlation (MSC);

Figure 12 is a graph of orientation error using frequency and model-based estimates;

Figure 13A is a graph of translation error with different pattern rotations using WPC;

Figure 13B is a graph of translation error with different pattern rotations using MSC; Figure 14A is a graph of tracking error for small translations using WPC;

Figure 14B is a graph of tracking error for small translations using MSC;

Figure 15 is a simulated grid with varying transmissions;

Figure 16A is a graph of the mean displacement error for focus patterns of varying transmissions using WPC; Figure 16B is a graph of translation error versus signal-to-noise ratio (SNR) for focus patterns of varying transmissions using WPC;

Figure 16C is a graph of mean displacement error for focus patterns of varying transmissions using MSC;

Figure 16D is a graph of translation error versus SNR using MSC; Figure 17 is an image of a simulated square pattern at different depths;

Figure 18A is a graph of translation error versus depth using WPC;

Figure 18B is a graph of translation error versus SNR using WPC;

Figure 18C is a graph of translation error versus depth using MSC;

Figure 18D is a graph of translation error versus SNR using MSC; Figure 19A is a graph of orientation error versus SNR for varying transmission,

Figure 19B is a graph of means error versus SNR for varying z position; Figure 2OA is a sequence of pattern images obtained at different illumination levels;

Figure 2OB is an image of segmentation of the image plane with pixels in medium gray and in bright gray; Figure 2OC is a graph of dark intensity versus mean bright intensity obtained from 200 images at different intensity;

Figure 21 is a series of images of a pattern being translated by a MCL nano-positioning stage;

Figure 22 is a series of images of a pattern cropped from a single high- resolution image used to simulate stage motion;

Figure 23 is a micrograph of an octopus muscle and a grid pattern;

Figure 24A is a series of experimentally acquired images;

Figure 24B is a series of simulated images created to match the experimentally acquired images of 24A; Figure 25 is an image mosaic assembled from 20 microscopy images;

Figure 26A is an image showing a frequency-focus slice of magnitudes;

Figure 26B is an image of the optical transfer function (OTF) as a height map with scaling;

Figure 27A-27D are images and graphs of focus versus frequency generated using a weighting function;

Figure 28A is an image showing a simulated sinusoidal pattern using four frequencies along each pattern axis;

Figure 28B is a graph of focal search;

Figure 29A is a graph of RMS tracking error as a function of pattern spacing;

Figure 29B is a graph zoomed in on a region marked in Figure 29A;

Figure 3OA is an image of a specimen with sparsely-distributed contrast;

Figure 3OB is an image of a tracking pattern used to track features of interest in the specimen of Figure 3OA; Figure 3OC is a graph of the SNR estimated for tracking the pattern in 3OB in the presence of specimen in 3OA;

Figure 31A is a graph of RMS tracking errors for tracking the optimized pattern at different distances from the in focus specimen; Figure 31 B is a graph of the same RMS tracking errors as a function of mean pattern SNR, considering the frequencies in Figure 3OB; Figure 32A is a micrograph of a section of frog brain tissue; Figure 32B is an image of the pattern used to track features of interest in the pattern of Figure 32A;;

Figure 32C is the SNR estimated for the pattern in Figure 32B in the presence of the specimen in 32A as a function of pattern distance from the focal plane;

Figure 33A is a graph of RMS tracking errors for tracking the optimized pattern in the presence of the tissue at different distances from the in focus specimen;

Figure 33B is a graph of the same RMS tracking errors as a function of mean pattern SNR, considering the frequencies in Figure 32A;

Figure 34A is a graph of mean focus estimation errors for simulated images of a sinusoidal pattern;

Figurer 34B is a graph of mean focus error versus SNR; Figure 35A is a graph the mean focus error versus z position; Figure 35B is a graph of mean focus error versus SNR; Figure 36A is a graph of the z tracking error for a non-optimized sinusoidal pattern in the presence of the frog brain section;

Figure 36B is a graph of mean focus error versus SNR; Figure 37A is a graph of focus error plotted against z position; Figure 37B is a graph of focus error plotted against SNR; Figure 38 is high-definition image mosaic of a large field of view.

DETAILED DESCRIPTION

I. Tracking with Structured Illumination Microscopy Structured illumination is a computer vision technique for obtaining information about an image scene. In this technique, a structured light pattern illuminates a scene in a controlled manner, and the effect of this illumination is used to determine scene properties, such as for example, geometry[1 ], [2], [3]. Novel applications of structured illumination to determine three-dimensional (3D) positions of a microscope stage are discussed herein. Figure 1 shows in one embodiment, a system for microscopy tracking, generally designated 100. System 100 can comprise a microscope 104 with a microscope stage 106. Microscope stage 106 is the platform on which a specimen sits during observation. Determining and tracking precise positions of stage 106 can be important in many microscopy applications. Stage tracking enables monitoring motion over large ranges, and uses a fixed stream of reference for moving or non-moving specimens. Microscope stage 106 can be repositioned to track the specimen, but the motion of stage 106 needs to be known independent of specimen motion to discern the complete motion of the specimen over an entire experiment. Given stage 106 position information allows measurements spanning multiple fields of view (FOV), for example to determine object size or measure the distance between feature landmarks. Tracking the X, Y, and Z position of stage 106 extends this capability to 3D. Stage tracking also provides information about the location in a specimen an observation is being made. A collection of smaller, high magnification images with accurate stage positioning data can enable assembly of a microscopy mosaic image that presents both fine structural details and broad context.

Using the structured illumination approach, a specimen holder or covering member 102 for holding or covering the specimen being examined under microscope 104 is used. Specimen holder or covering member 102 can comprise for example, a slide or slide cover and can be located above or below specimen. A pattern can be formed in or on specimen holder or covering member 102 to modulate the light transmitted through, for example, a bright field microscope. The pattern acts as a transmission mask, partially attenuating the light illuminating specimen in a structured manner. An imaging device 108, such as a digital camera, can collect images from microscope 104. Analysis of pattern appearance by using an image analyzer 110, which may be implemented in software, can enable tracking the stage in 3D, that is laterally in X and Y directions and axially in Z direction. For example, phase-based tracking on the frequencies present in the pattern can enable lateral tracking. Defocus analysis on the pattern frequencies can enable axial tracking. This method can enable three-dimensional stage tracking independent of specimen motion. Image analyzer 110 can also be equipped to remove pattern from the image once an image is captured. As a result, pattern free images may be used to build high-resolution, wide field of view microscopy mosaics.

Figure 2 illustrates exemplary overall steps of a method for microscopy tracking, generally designated 200, according to an embodiment of the subject matter described herein. Referring to Figure 2, in step 202, the method can include providing a specimen holder or covering member for a microscope where the specimen holder or covering member may comprise a slide or slide cover. The specimen holder or covering member can be positioned above or below a specimen to be viewed using the microscope. In step 204, images of the specimen with the patterned specimen holder or covering member can then be captured using the microscope. In step 206, the images are analyzed, and location of the stage is determined through at least one of frequency, phase, and defocus analysis of the pattern. In one example, tracking stage position can include determining absolute stage position using frequency and phase information encoded in the structures in the pattern. For example, different lateral axes in the pattern may have structures with different spatial frequencies. The phase of the pattern along each axis in captured images may provide lateral stage coordinates along each axis. Each position in the pattern could have a unique combination of phases, so this coordinate system could indicate absolute stage position.

Tracking of locations of features in the lateral directions includes tracking the locations using spatial frequency data from the pattern. Tracking in the axial direction includes using defocus analysis. This can determine position from, for example, a focal plane. The tracking can be performed by using image analyzer 110, which may be implemented using a computer having imaging software using the algorithms described below. In step 208, pattern removal may be performed by image analyzer 110, for example, by using various forms of occlusion removal as described below. In on example, in step 210, a mosaic of images may be created using the pattern. For example, plural images obtained at a first field of view may be assembled into a mosaic image of a second field of view using the pattern. Examples of such mosaicing of images will be described in detail below. Aside from mosaicing, other applications for this novel system can include gaining absolute stage position data from encoded frequency and/or phase data, measuring samples across several fields of view, and analyzing live specimens as stage movement can be independent from specimen movement.

II. Theory of Stage Tracking

In a bright-field microscope, light from a uniform light source is transmitted through a specimen and collected by a series of lenses to form an image on a camera sensor. Image contrast is generated by the transmission and diffraction of light through the specimen. An image of the specimen can be predicted by geometrical optics and modeled as a thin, planar, transmission mask S(x; y), the image formed by the microscope is: l(x,y) = L(x,y)[h(x,y) ® S(x,y)] where L(x; y) is the illumination at the focal plane, h(x; y) is the microscope's point-spread function (PSF) at a focal plane, and ® S(x,y) is the convolution operator.

A small pattern, P(x; y) can be inserted into the microscope's light path on the specimen holder or covering member, slightly above the specimen. The pattern thus acts as a transmission mask that attenuates the illumination field, but its image is convolved with a different part of the microscope's 3D PSF. The image formed on the imaging device becomes the multiplication of images of each plane according to the following approximation:

I(x, y) = L(x, y)[h(x, y, Z₈ ) ® S(x, y)] [h(x, y, z _p ) ® P(x, y)] where z_s and z are axial positions of the specimen and pattern respectively.

Note that this model does not precisely account for the diffraction of light between the pattern and specimen. However, because the distance separating these planes is much smaller, for example as little as 5 μm compared to the distance the light diffracts within the microscope, for example 160 cm, this approximation is reasonable. This approximation gains the additional benefit that the ordering of the pattern and specimen layers is unimportant. The position of the pattern is fixed with respect to the slide, which remains fixed with respect to the microscope stage. The pattern acts as a semitransparent transmission mask, and it moves with the microscope stage. Introducing the pattern, then, transforms the stage tracking problem into a pattern tracking problem. Obtaining accurate stage tracking requires that the pattern be tracked independently from the specimen, and a successful tracking method will assume little about the specimen while obtaining accurate pattern tracking. Specifically, pattern tracking holds an advantage over any tracking method that relies on specimen information because it applies whether or not the specimen has sufficient contrast and whether or not the specimen is moving. Because the pattern can be semitransparent and has a known appearance, the effect of the pattern can thus be removed using an extension of the stationary occlusion removal techniques described further below.

The pattern tracking problem can be broken down broadly to a method of finding the orientation and size of the pattern in all input images, determining the axial position (focus) of the pattern with respect to the specimen, determining the lateral position of the pattern along each pattern axis independently or simultaneously, and repeating this iteration for any additional information gained at each step.

III. Determining Lateral (X, Y) Pattern Orientation

Figure 3A depicts an image of a regular pattern and Figure 3B depicts its Fourier transform. The regular pattern can be composed of two identical, orthogonal square patterns. Note that this pattern's Fourier transform has a fundamental frequency determined by the spacing between square centers plus harmonics along each of the pattern's coordinate axes. The harmonics occur at integer multiples of the fundamental frequency. The first step in stage tracking is to determine the orientation of the pattern's coordinate axes so that tracking can proceed along each axis independently.

An initial estimate of the pattern orientation is determined with a search in frequency space. Rotation in the spatial domain corresponds to rotation in the frequency domain, so the frequency components from the pattern fundamental and harmonics lie a set distance from the center of frequency space. Figure 4A illustrates a subset of Fourier transform magnitude for the pattern in Figure 4B. A search for maximum magnitude along frequency bands of the fundamental and harmonics marked as rings in Figure 4A for the square wave pattern in Figure 4B yields four peaks for each frequency, spaced 90° apart. The orientation estimate can then be refined through a one- dimensional (1 D) optimization along the frequency band. The mean of all peak orientations for all target frequencies provides an orientation estimate. This search process can be adapted for any other pattern, given the frequencies and orientations present in the pattern.

If pattern spacing is only roughly known, an orthogonal search provides a refined estimate for pattern spacing given an initial orientation estimate. In this case, the 1 D optimization operates along the pattern axis to find the frequency with maximum magnitude. Figures 4C and 4D show the magnitude along each axis U and V respectively, with vertical lines marking maxima. The orientation and frequency searches alternate until convergence.

IV. Pattern Tracking Because the specimen holder or covering member orientation stays fixed during stage motion, one can use multiple images to estimate pattern orientation. Once orientations are known, tracking requires recovering the translation of pattern from frame to frame. Several image analysis techniques can recover translation from images, such choices for pattern tracking include feature tracking, spatial correlation-based matching, and phase correlation [4]. The fundamental image analysis task in tracking patterns is to separate image information in pattern layers from specimen layers. Spatial correlation and phase correlation are discussed in detail below, as they are not as susceptible to interference from specimen layers as feature tracking.

Model-based Spatial Correlation (MSC)

Image registration is a class of image analysis techniques for determining motion from images. Assume two (nxm) images, I1(x, y) and I2(x, y), are taken of specimen oriented perpendicular to the optical axis. If there is only translation of the object (or the camera) perpendicular to the optical axis between the capture time of the two images, the images are related by I1(x, y) = I2(x - u, y - v). This suggests that finding the translation (u, v) that describes the translation between the two images amounts to minimizing an energy function, such as the L2 distance:

E(U₁ V) = ^ [I₁(X₅ Y) - I₂(X - U₅Y - V)]² χ»y

An image comparison metric known as normalized cross-correlation (NCC) can be derived, for example from the expansion of the image comparison term in the equation above [5]. For two image patches I _x (x, y) and 1₂ (x, y):

d² = E (I₁(X₅Y)² + I₂(X₅ y)² - 2I₁(X₅Y)I₂(X₅Y)). χ,y

Assuming the terms El_j(x,y) and El₂(x,y) as approximately constant χ,y χ,y and normalizing the remaining image comparison term yields:

where T(x, y) = l(x, y) ∑I(^χ ₅y) . ^ar)d ^a" sums are over the image plane. nm Because NCC normalizes the distribution of intensities in the compared images, it compares the landscape of the two images and is robust to intensity shifts. Note that finding the translation between two images requires finding the alignment that maximizes the NCC metric equation above.

Image registration provides another method for tracking the pattern layer in structured illumination microscopy. A model of the pattern's appearance can be used to generate simulated pattern images. An optimization procedure finds model parameters to maximize the NCC metric between a microscopy image and the simulated pattern images. That is, given a model that produces a pattern image, M(u, v)-> Pu,v(x, y), find translation parameters (u, v) that minimize

E(u, v) = - NCC(I(X, y); Pu,v(x, y)). Hence, this method will be referred to as model-based spatial correlation (MSC).

The image formation model needs to account for pattern design, orientation, translation, and the effect of defocus. For a fixed orientation and focal depth, finding the translation of pattern in an image involves a search over the parameter space of all possible translations. An exhaustive search that considers sub-pixel displacements is slow, so a practical approach relies on maximizing the NCC objective function using an optimization strategy. The shape of the objective function depends on the pattern and the specimen.

In one embodiment, the optimization method can involve a coarse- resolution search through parameter space to initialize a Nelder-Mead simplex- based optimization [6]. For 2D optimization, this technique maintains a set of three connected points in parameter space-the simplex-and the values of the objective function at those points. The simplex is updated by evaluating the objective function near the simplex, moving the high-valued simplex vertex to lower valued locations in an attempt to scale towards the minimum. This approach, of course, depends on having a good initial guess to initialize the simplex. The initial search proceeds by first selecting an initial parameter range to allow for the largest translation of pattern between frames that can be distinguished. Then a target microscopy image can be smoothed with, for example, a Gaussian filter and down-sampled to a coarse resolution to make the image simulation and NCC comparisons fast. A down-sampling factor that is a power of two is chosen to ensure that a single pattern repetition covers at least an 8x8 pixel region in the down-sampled image. The standard deviation of the Gaussian filter is chosen to be the square root of the down-sampling factor to avoid aliasing in the down-sampled image. An initial translation step size can then be chosen to ensure that each pixel within a down- sampled pattern region is sampled at least once. Coarse pattern images can be simulated for each translation using large virtual CCD sensors to account for the reduced resolution. A NCC metric can be computed for each translation to assemble a correlation map over the parameter space. The parameter range and step size can be contracted around the maximum in the correlation map by a factor of two. The former procedures can be repeated twice at the same resolution and using a higher-resolution target image and smaller CCD sensors in the simulation. When the search has completed on an image resolution half the size of the original image, the current translation and range initialize a modified Nelder-Mead simplex-based optimization. The Nelder-Mead optimization procedure usually conducts an unbounded search in the parameter space, but setting the objective function to a very poor value outside the range determined in the multiple resolution searches constrains the parameter space. The initialization method arrives at a translation that is no more than 1 pixel away from the final optimized translation. This spatial correlation method is evaluated for structured illumination microscopy patterns.

A similar model-based correlation approach can improve the orientation estimate provided by frequency analysis. Given an initial orientation and position estimate, pattern images can be simulated at different orientations to find one that maximizes the NCC metric. The 1 D parameter search with an initial estimate close to the optimum is readily implemented, for example, with a golden section optimization [7].

Phase Correlation Phase correlation can provide an efficient registration method for recovering translation transforms, such as those introduced by moving a microscope stage [8]. Phase correlation depends on the Fourier shift property, for example, a translation in the data domain is equivalent to a phase shift in the frequency domain. In one embodiment, images 11 and I2, for example: I2(x; y) = M(X - xO, y - yθ), the cross power spectrum yields:

(U_{5 V}) _e-i(uxo+vyo)

where l(u, v) is the Fourier transform of l(x, y) and I^* is the complex conjugate of I. The inverse of the cross power spectrum is a Dirac delta function centered at translation coordinates (xO, yθ) for example, which can be found by locating the maximum in C(x, y), the Fourier transform of C(u, v). Efficiency is one of the major advantages to phase correlation. The method requires three (3) Fourier transforms, two pixel-wise image products, a pixel-wise image division, and a search for a maximum. The computational complexity of phase correlation is therefore limited by the Fourier transform, which for a pair of images with N pixels has O(N log N) complexity [9]. In contrast, spatial correlation has O(N²) complexity at potentially every translation, a pixel-wise image product and summation performed.

There are two barriers to using phase correlation to track a regular pattern. First, because phase correlation treats all frequencies as equal contributors to the registration, this method can require broad spectrum support to recover a translation [10]. Second, phase correlation requires applying a windowing function to the image data to attenuate discontinuities at image edges. This weights information in the center of the image more heavily than information at the edges of the image. Though this is less of a problem when considering rigid translations of objects covering the whole image plane, in pattern analysis it is desirable to treat all valid pattern data as equal contributors to registration. Specifically, it is unknown where in the image the pattern information might dominate the specimen information; if this occurs only near the borders of the image for example, windowing would discount this important tracking data.

The phase estimate for particular pattern frequency should be very accurate to obtain a reasonable tracking accuracy. For example, a phase error of Eφ measured in degrees, the pixel tracking error is

Epx = Eφ

36Of where f is frequency of the component used for tracking. For a fixed phase error, tracking accuracy can improve for example, with higher frequency and shorter wavelength pattern components. The range of displacements between frames that can be unambiguously determined with a single component, however, can be limited to half a wavelength in each direction. The following modifications to phase-based tracking are useful for example, in specific imaging applications discussed wherein tracking may involve only the pattern. A review of the Fourier transform and discrete Fourier transform (DFT) provides context for discussing attempted approaches to overcome potential barriers to pattern tracking with phase correlation.

Discrete Fourier Transform (DFT) Fourier transform is a representation of a signal as a collection of frequency components. For a continuous, infinite 1 D signal, f(x), the component of the Fourier transform for a frequency u is given by: oo

F(u) = J f(X)₆-¹²⁷¹¹¹Mx₁

—oo where i = v—1 ■ The Fourier transform is for example, a change in data representation between the spatial domain (all spatial positions x) and the frequency domain (all frequencies u). Each component of the Fourier transform for example, has a magnitude that can indicate energy of the frequency component contained in the signal and a phase that indicates the relative position of the frequency component in the signal. For a complex component, u, the magnitude and phase are

M(u) = |F(u)| = V5H(f(u))² + 3(F(u))² ,

where in one embodiment, SR and 3 represent the real and imaginary components of a complex value, respectively. The magnitude of a signal is often reported in decibels, a unit-less measure that compares for example, component magnitude M to some reference magnitude M₀ on a logarithmic scale. For images, the scale convention is

The DFT is for example, a discrete approximation to the Fourier transform. For aa 11DD ddiissccrreettee ssiiggnnaall,, ff((xx)),, ddeeffiinneedd over the domain x = [0...N - 1], the k* component of the DFT is given by N-I -j^

F(k) = ∑ f(x)e ^N k=0...N - 1. x=0

Each component of the DFT is the correlation of the input signal with a sinusoidal basis function (from Euler's formula, e^"1 = cos θ - i sin θ ). As the basis functions are periodic and have infinite extent, the DFT treats the input signal periodic as well. As such, the DFT of a signal computed over an infinite signal domain is equivalent to the DFT of an infinite number of copies of the signal computed over an infinite domain. For example, consider the finite signal as a truncated portion of an infinite signal. The finite signal can be created by multiplying the infinite signal by a rectangular windowing function that has value one within a finite domain and zero everywhere else. In fact, there are many forms of windowing functions used in Fourier analysis. Generally defined, a windowing function has value zero everywhere outside a finite domain. It can be said for example, that every DFT computation on measured data involves windowing because every observable signal is finite in extent.

From the convolution theorem for example, multiplication in the spatial domain is said to be equivalent to convolution in the frequency domain [11]. The effect of multiplying a signal by a windowing function in computing the DFT is convolution in the frequency domain by the Fourier transform of the windowing function. The result of this convolution is that a single frequency component can include contributions from neighboring frequencies, an effect often known as leakage. Selecting a windowing function can involve a trade off between the width and fall-off of the windowing function's Fourier transform [12]. These parameters can affect how many neighboring frequencies contribute to a single component and the magnitude of each contribution, respectively. In one embodiment, an ideal windowing function would have a frequency domain representation of a single point, but this would require a windowing function of infinite extent in the spatial domain.

Harmonic analysis-determining magnitude and phase of the frequencies in a signal-is an estimation problem [12]. Figure 5 illustrates the effect of windowing on parameter estimation of sinusoidal signals of the form f(x) = Acos(2 π ux + θ )

In Figure 5A, A = 10, u = — , and θ = 90°. This signal for example, is

N composed of a single component that has an integral number of periods over the data domain x = [0...N]. The dotted lines in the signal plot show how the signal repeats outside the window extent, as considered by the DFT. The magnitude is then plotted in decibels, with the reference magnitude for all plots in Figure 5 the maximum magnitude from the signal in Figure 5A. Because

]£ there is only a single frequency (u = — ) in the signal and all of its energy falls

N into a single DFT bin, the magnitude plot for this signal exhibits two sharp k spikes corresponding to frequencies ± — . The magnitude and phase

N computed for this DFT component are A = 10.0 and θ =90°, matching the input signal precisely.

Figure 5B presents the result obtained when a window is applied to the signal from Figure 5A. Here for example, a Blackman-Harris window is used because it provides a good balance between the width and fall-off of the frequency domain representation [12]. The effect of spectral leakage though is evident in this case-energy from the single frequency component in the signal spreads into neighboring frequencies. The signal parameters estimated in this case are A = 10.0 and A = 90°, matching the input signal precisely. Figure 5C for example, presents a signal composed of a single component that has a non-integral number of periods over the data domain-the

5.125 function has the same form as that of Figure 5A with u = . The dotted

lines show how the signal repeats outside the window extent, as considered by the DFT. Although the frequency present in this signal differs only slightly from that of Figure 5A, the non-periodicity at the data boundaries introduces a step edge to the signal, which has energy at all frequencies, as seen in the magnitude plot. The magnitude plot shows that the peak amplitude is decreased slightly (A = 9.64), but energy from the single frequency component present in the signal leaks energy into several neighboring DFT bins. A consequence of phase correlation, however, is that the phase shift measured at the peak amplitude can be grossly mis-estimated at θ = 112.7°.

Figure 5D presents the result obtained by applying a Blackman-Harris window to the signal from Figure 5C. In the frequency domain, the windowing function attenuates the signal to be continuous at the boundaries, eliminating the step edge seen in Figure 5C. As such, peak magnitude and phase are estimated to be A = 9.94 and θ = 112.5°, respectively. Although the windowing has removed the discontinuity at the data edges, the phase estimate can remain inaccurate as input signal does not exactly match one of the DFT bins.

Weighted Phase Correlation (WPC)

The use of different windowing functions can improve determination of whether a signal is present [12]. However, windowing may not necessarily improve the phase estimation task. Phase estimation is crucial in pattern analysis, therefore, frequency-specific phase analysis technique weighted phase correlation (WPC), is presented. WPC computes individual components of the Fourier transform for the frequencies present in a pattern.

Recall that the discrete Fourier transform (DFT) is the correlation of a signal having a set of sinusoidal basis functions. The DFT for a 1 D signal l(x) can be represented by:

N-I j^2πkx l(k) = ∑ I(x)e ^N k=0... N-1 , x=0

Target frequency fu can be defined by wave number u such that fu = — . Also,

N

assume wavelength, λu = — , can be an integer (λu Θ Z). As l(x) can be a

signal containing only frequency fu, then DFT l(k) can have single magnitude spiking at k = u if f _u cycles an integral number of times in the data domain.

Alternatively, the phase at frequency fu in the signal can be computed through a reweighting of the input signal, which is a goal of WPC.

As l(x) is periodic at intervals of λu, the partial summation: ^J λu

P λ(x) = £ P(x + jλu) , x = 0...λu - 1 , j=o can align components of l(x) that come from the same position along the wavelength and would be multiplied by the same part of the Fourier basis function in l(k) immediately above. The number of terms from the input signal that contribute to each point on Pλ (x) can be described by:

λu

Wλ (x) = ∑ l , x = 0...λu - 1. j=o

A reweighted signal,

Vλ(x) is an appropriately scaled representation of l(x) over a single cycle. At this point, a windowing function can be applied to the reweighted signal. The single frequency component at the target frequency fu can be computed as λu-l -£≡

Y(U) = ∑ Iλ(x)e ^N _t x=0

The phase component of l'(u) is for example, the same as the phase of fu in the input signal. The magnitude component differs by a scaling factor dependent on the number of repetitions of the signal over the original data domain. Figure 6 depicts WPC for a 1 D signal with a single frequency

4.5 component fu = — . For example, 6A is the original input signal; 6B is the

N partial sum of 6A by target wavelength. Figure 6C weights showing the number of contributions for each element of Figure 6B, and Figure 6D is the partial sum of Figure 6B reweighted by Figure 6C. Figure 6E is the Fourier basis function with target wavelength, and Figure 6F is weighted signal 6D multiplied by

Fourier basis function 6E.

This formulation has assumed that the wavelength of the target frequency is integral so that the correlations can proceed as component-wise multiplication and summation. If for example, λu is non-integral, the input data can be resampled to enforce this condition. The proper resampling technique chooses a resampling interval and determines precisely how much energy from each initial signal sample falls within each resampling bin. A simple resampling approach chooses a resampling interval that is slightly larger than the original data interval so that each input sample point maps into two resampling bins, followed by reconstruction of the signal data with a triangular kernel (linear interpolation).

Extending WPC to 2D requires accounting for the orientation of the u v target frequencies in the image, (fu = — ; fv = — ). For an axis-aligned

N N pattern, each axis can be treated independently. For example, v = 0 for u ≠ 0 and u = 0 for v ≠ 0. An example of 2D WPC appears in Figure 7. The first step in 2D WPC is to axis-align the pattern by rotating the image as in Figure 7A. The rotated image contains regions of invalid pixels, where the image data is unknown; these pixels are discounted from the weighting function (these regions are visible, for example at the top and bottom of Figure 7D). Any pixels that receive no contributions during the partial summation are discounted from the reweighting step because reweighting is undefined at those locations. Otherwise, the 2D computation follows very closely from the 1 D computation. Figure 7A is the axis-aligned input image with target wavelength marked by dashed vertical lines. Figure 7B shows partial sum of 7A by target wavelength.

Figure 7C shows the number of contributions for each element of 7B. Figure

7D is the partial sum of Figure 7B reweighted by Figure 7C. Figure 7E is the

Fourier basis function with target wavelength. Figure 7F is the weighted signal of Figure 7D multiplied by the Fourier basis function of Figure 7E. Note, however, that the reweighting in WPC may not avoid windowing.

Although it addresses the problem of the target frequency not falling within a single DFT bin, it does not for example, address the problem of other frequency components in the input signal that repeat a non-integral number of times in the windowed domain. The data has been aligned and reweighted according to a particular frequency that is present in the signal, but there are many other frequencies present which can still produce leakage in the frequency domain. A windowing function can be employed as usual to attenuate the effect of discontinuities at the data boundaries. However, any fractional wavelength segment of the input signal may introduce discontinuities in the middle of the data domain. This can be particularly problematic for rotated images, which have discontinuities along the diagonal boundary between the signal domain and the unknown data domain.

Comparing the phases for pattern frequencies in two images yields a phase shift corresponding to the translation of the pattern in the images. The phase shift can be measured for any frequency present in the pattern along each axis. Phase shifts from multiple frequencies can be combined to improve the translation estimate, with proper consideration that the phase shifts from higher frequencies may indicate multiple cycles over the translation.

Computing the translation between a set of images and a reference image can provide translations that register all images to a common reference frame.

In one embodiment, the process described can be thought of as similar to the computation of particular DFT components. However, the computation is carried out for a specific set of frequencies with potentially non-integral wave number, and a reweighting is applied to account for this and for missing data in the input image. The result of this discrete, frequency-specific, weighted correlation can be used to compute the phase for any frequencies that provide information about patterns in an input image.

Comparing DFT and WPC in 1D

Estimating phase parameters using general phase correlation is not usually a problem for wide bandwidth signals, such as for example when a complete image undergoes translation. In this case, phase correlation relies on corroborating measures from many frequencies in the full DFT to arrive at an average phase shift. For narrow bandwidth applications however, the phase shift estimation is performed on a small number of samples in frequency space, and for each sample there is significant prior knowledge about the signal, for example, its periodicity is known. The WPC method relies on this prior knowledge and removes the requirement that the signal has an integral number of repeats in the image. Comparative evaluation of the phase estimation task for 1 D signals of different frequencies using DFT and WPC provides further justification for the WPC approach. Following an algorithm evaluation approach outlined by Trucco and Verri [13] a set of signals representing a range of test scenarios are generated. Each signal contains a single frequency component with a known amplitude and phase shift in the form: f(x) = Acos(2 π ux + θ )

For each test scenario for example, 20 instances of the signal are corrupted with zero-mean, Gaussian-distributed noise with standard deviation 10% of the signal amplitude. The phase shift is estimated using both the DFT and WPC techniques, and the mean phase estimate can be computed from the 20 samples for each test scenario. The signal frequency varies over the range {u lλ u = — , λ = [5, 6,...N]} and phase varies over θ = [-π ... π]. In one

embodiment, the results, the root mean squared (RMS) error and standard error are computed to summarize the behavior of each method. The RMS error is represented by:

Ei = xi - x_τi

where xi are the N computed mean parameter estimates and xTi are the known ground truth values. The standard error is:

where E can represent the mean of all computed errors [7]. The RMS error and standard error can thus provide an estimate of accuracy and error in the accuracy estimate, respectively. Figures 8A and 8B plot the phase estimation error for each estimation method using a rectangular window. Phase estimation for a target wavelength using DFT and WPC methods on a signal containing one frequency is shown. WPC performs well throughout the range of tested frequencies, u = [5...N]. DFT, however, only performs well at frequencies that are periodic in the window domain, {u u = — , k = [1... N - 1]}. Figure 8A shows the mean phase error

N for all methods using a rectangular window, DFT with the nearest neighbor and linear interpolation. Figure 8B shows the mean phase error from 8A for the WPC method only. Note the change in the y axis scale.

Table 1.0 is a comparison of phase estimation using DFT and WPC on 1 D single frequency signals with 10% Gaussian-distributed noise. Nearest neighbor and linear interpolation are used to estimate frequency components that fall between DFT bins, as indicated.

Table 1.0 records the RMS error and standard error for phase estimation using DFT and WPC with several popular windowing functions. Parameters are provided for some windowing functions, using the terms provided by the Harris [12]. These parameters adjust the width and fall off of the functions' Fourier transforms. In each case, the best of several parameter trials are reported. None of the DFT methods reliably estimates phase parameters over the entire frequency range, regardless of windowing function. WPC, however, performs considerably better, with the best result obtained with the rectangular windowing function (no windowing). This evaluation highlights one of the shortcomings of using the DFT for phase estimation using a small set of frequencies. For example, there are only a few frequencies for which the estimation will be accurate. The sections that follow discuss characteristics of well-designed patterns for structured illumination microscopy, but to foreshadow a little, it is desirable to have a number of long wavelength components in the pattern, as this enables tracking over longer ranges and with greater degrees of defocus. WPC enables the use of more wavelengths for which phase estimation will be accurate.

The situation changes, however, in the presence of other frequency components. Figures 9A and B show the phase estimation error for a signal of the form: f(x) = Acos(2 π ux + θ ) + 0.1 Acos(2 π ucx).

That is, the single frequency signal computed from above has a corrupting signal added with 10% amplitude. In this case, the corrupting signal's wavelength is λ = 56.5, and therefore does not align with the target frequency in any test case. Figure 9A is the mean phase error for all methods using a rectangular window, DFT with nearest neighbor and linear interpolation are nearly identical. Figure 9B is the mean phase error from Figure 9A, for the WPC method only. Note the change in y axis scale. Figures 9 illustrate the ability of WPC to estimate phases at most any frequency and are degraded by the presence of this corrupting signal, and is significantly more degraded at frequencies near the corrupting signal's frequency. The maximum error in this situation is approximately 6°, or 1.7%. For example, Table 1.1 is a comparison of phase estimation using DFT and WPC on 1 D signals containing two frequencies with 10% Gaussian-distributed noise. Nearest neighbor and linear interpolation are used to estimate frequency components that fall between DFT bins, as indicated.

Table 1.1 records the RMS error and standard error for phase estimation of the corrupted signals using DFT and WPC and the windowing functions from Table 1.0. The results illustrate how windowing does not improve phase estimates of single frequency components. While the accuracy of WPC may not be as good in the presence of corrupting signals, in conclusion, it performs better than DFT. Whether this is sufficient for tracking patterns in structured illumination microscopy depends to a large extent on the specimen being observed as discussed later.

Evaluation Using Simulated Images

A combination of simulation and real experiments are outlined which evaluate the two lateral microscope stage tracking approaches of weighted phase correlation (WPC) and model-based spatial correlation (MSC). The experimental data uses a square grid pattern printed on a cover slip. Simulated images of this pattern undergoing translation are used to evaluate the tracking methods using known ground truth displacements. The simulated pattern is based on a square transmission electron microscope (TEM) grid composed of a set of orthogonal 6 μm bars spaced 25 μm apart. This arrangement results in a set of regular square regions with a 19 μm hole in the center. The simulation is designed to mimic the way the pattern used in real experiments is imaged in the microscope, considering such parameters including the following: β The physical size of the grid pattern (25 μm) β The objective lens used for observation (4OX, 0:65NA) * The optical transmission of the pattern β The orientation of the pattern β The axial position of the pattern β The size of the CCD sensor elements (9 x 9 μm) » The shot noise from image acquisition (Poisson noise with 100 photoelectron per count ADC gain)

Figures 10A-C show the steps used to simulate imaging of a structured pattern with transmission α=0.5. In one embodiment, the first step in simulating an image of the grid pattern can be computing the image formed on a virtual microscope's CCD predicted by geometrical optics. The pattern is given a transmission α and is illuminated by a uniform light source of intensity L. An image of the pattern is generated by evaluating an analytical function of the grid

with a specific orientation θ and translation [dx, dy] with 4 x 4 oversampling for each CCD element. An example of the image predicted by geometrical optics is shown in 10A. This supersampled image is convolved with the point spread function predicted by Fourier optics for the grid displaced from the microscope's focal plane by a distance dz. An example of the image predicted by Fourier optics is shown in Figure 10B.

The Fourier optics image is down sampled to obtain the ideal intensities read by each CCD pixel. Noise is added to each pixel to model shot noise using either a Gaussian or Poisson noise distribution. The Poisson noise model simulates the shot noise seen in a real camera. The number of photons arriving at a sensor site is governed by a Poisson random process [15]. In image simulation, a parameter modeling analog-to-digital converter (ADC) gain (photoelectrons per count) converts ideal intensity values to photoelectron counts at each pixel. A Poisson-random number with variance equal to the electron count is generated at each pixel and converted back to an intensity value. With high ADC gains (high photoelectron/count ratios), however, the Poisson distribution can be closely approximated by a Gaussian distribution [14]. For the Gaussian noise model, Gaussian-distributed noise is added to each pixel with a standard deviation that is a percentage of the brightest image value. The noisy image is finally quantized to whole-integer intensity values to produce the final simulated image of the pattern, as seen in Figure 10C.

This process enables evaluating the tracking algorithm for varying pattern transmissions, rotations, translations, focus positions, lens parameters, and noise models. The following approach makes use of the approach outlined in [13]. For each test case for example, 5 images are generated with noise to simulate repeated measures. For each simulated image, the rotation of the pattern and the displacement from a reference image is measured using WPC and model-based spatial correlation (MSC). The difference between the mean parameter estimate and the ground truth (the mean error) indicates how accurately each tracking algorithm handles a particular test scenario. The standard deviation of measurements indicates the algorithm's robustness to noise.

Shot Noise

Camera shot noise was evaluated on the tracking algorithm. For example, in one embodiment the simulated grid can be given a transmission of α = 0.5, placed unrotated and in focus, and then translated between -25 and 25 pixels along the x-axis. The grid image can be corrupted with either Gaussian-distributed noise or Poisson-distributed noise as seen in Figures 11 A and 11 B. Figures 11 A and 11 B show the mean displacement error for different noise models with the two tracking methods. The left hand side charts for Figures 11A and 11B show Gaussian noise with standard deviation as a percentage of maximum image intensity. The right hand side charts for Figures 11 A and 11 B show Poisson noise with different ADC gains (electrons/count). Figure 11A is the error using WPC tracking, and 11 B is the error using MSC tracking. Note the change in y axis scale. The tracking accuracy for WPC is consistently high, independent of noise model, with higher levels of noise disrupting the accuracy more. The WPC method attains far more accuracy than the MSC method in this scenario. The MSC method is nearly invariant to noise levels in this range. These results also indicate that tracking error does not depend on translation distance.

Table 1.3 below illustrates RMS displacement error and standard error ( σ E) for different noise models using the two lateral tracking methods. Results are for tracking an in focus square grid pattern with transmission α = 0.5, and σs denotes the standard deviation of the shot noise.

Table 1.3 reports shot noise standard deviation, σs, for each scenario. For example, for the WPC model, similar shot noise levels produce comparable tracking errors using both the Gaussian and Poisson models. This indicates that the distributions are interchangeable for modeling shot noise, especially when high ADC gains are used. The RMS displacement error remains bounded by 0.010 pixel for all noise models using WPC, which corresponds to an error of 2.3 nm for the simulated microscope's configuration. Some caution should be used in interpreting this value, as this error may represent the scenario of imaging an in focus pattern with precisely known spacing and orientation and no corrupting specimen. As seen in the table and Figures 11A and 11 B, the MSC method appears worse in this comparison, maintaining a nearly-constant error of 0.25 pixel, or 56.3 nm. Pattern selection may improve this level of error as discussed later. The remainder of the system parameters fall into two categories: those that do not degrade the pattern signal (rotation and translation) and those that do (transmission and focus, see discussions to follow.

Rotation Pattern and camera sensors can be placed at any arbitrary orientation.

The tracking algorithm relies on being able to both accurately determine the orientation of the patterned specimen holder or covering member and track the specimen holder or covering member at any orientation. These tasks are evaluated independently, with two methods used for each. Orientation can be evaluated with the frequency-based estimation and the model-based refinement associated with MSC tracking. The model-based method can depend on initial estimates of both the pattern position and orientation while the frequency-based method does not. In this evaluation, the estimate provided by the frequency-based method is used to initialize the model-based method to reflect how the two methods are used in practice. The position estimate is provided by the known ground truth position of the pattern.

Figure 12 shows the mean orientation error in response to different pattern rotations obtained by the two methods for pattern rotations within θ = [0...90] ° for an in focus square grid pattern of transmission α=0.5. Initial orientation for the model-based estimation is provided by the frequency-based estimation. Error bars are too small to be resolved for the data in Figure 12. Orientations outside this range repeat for this square pattern. The data shows that using the frequency-based estimation, the orientation error is bounded within ±0.4° of the true value. The offset of a pixel at the edge of the 256 pixel test images whose orientation has been miscalculated by this amount is 128 tan(0.4) « 0.89 pixel (0.20 μm). Over all estimates, the RMS error using frequency-based estimation is 0:24° with standard error 0:073°. The model- based estimation provides nearly a factor of four improvement, with an RMS error of 0.066° and standard error 0:019°. This presents a strong argument for refining the orientation estimate after obtaining a position estimate.

Figures 13A and 13B show translation error obtained for tracking a grid pattern displaced with different rotations. Figure 13A illustrates the WPC translation error in x and y, while Figure 13B illustrates the MSC translation error in x and y. The mean translation errors obtained using the two tracking methods where the pattern orientations are provided from ground truth values and the pattern translation is dx = [10, 0]^τ pixel. The data reveals MSC outperforms WPC for most scenarios, with exceptions at 0 and 45°. The large errors with MSC are again attributable to the particular pattern being tracked. Table 1.4 summarizes these results along with tracking errors for small pattern displacements, discussed next.

Translation The previous test case demonstrated that the tracking method maintains accuracy at arbitrary orientations-even though the method is also dependent on estimating those orientations precisely. In an evaluation of the range of pattern translations, then, it suffices to investigate the translations in one direction for an axis-aligned pattern. The maximum frame-to-frame displacement along a pattern axis is limited to half the wavelength of the pattern's fundamental frequency along that axis. A greater displacement will alias to a displacement in the opposite direction, so larger displacements cannot be recovered accurately without additional knowledge. The error for large displacements of the pattern was already investigated during the discussion of shot noise. The RMS translation errors for the Poisson-distributed noise with ADC gain of 100 photoelectrons / count are reported in Table 1.4 below. Table 5.4 illustrates RMS error and standard error in translation estimation using WPC and MSC for different pattern orientations and small translations. Gross errors have been discounted from the MSC orientation estimates. All measurements are given in pixels.

Because the WPC tracking algorithm is based on phase comparisons and MSC is based on matching simulated images to an observed image, there is nothing that locks either algorithm to measuring whole-pixel displacements. Figure 14A and 14B show translation error using WPC and MSC in each direction for % pixel displacements along the x-axis in the range dx = [-2...2] pixel. Figure 14A uses WPC and 14B uses MSC. Vertical error bars indicate the standard deviation for each measurement taken over 5 samples. Table 1.4 records the RMS errors for these test cases. The data reveals that WPC vastly outperforms MSC, maintaining a worst-case error of 0.025 pixel (5.6 nm). The optimization of MSC, on the other hand, appears to be trapped in a local minimum, always reporting zero translation, which gives rise to the obvious structured error in Figure 14B. Once again, this is a problem that the square grid pattern poses to the optimization strategy; hence pattern selection will be discussed below.

Transmission

Figure 15 shows a simulated grid pattern with varying transmissions, from left to right, α= {0.45, 0.60, 0.70, 0.85, and 0.95}. The remaining test parameters degrade the quality of the signal used for tracking. The pattern transmission and lamp intensity together determine how much information is available in an image sequence for phase-based tracking. Optical transmission is the percentage of light that passes through a semitransparent object. The grid pattern acts as a transmission mask, segmenting the image plane into bright and dark regions. Figure 15 shows simulated images of square grid patterns with varying transmissions and constant lamp intensity.

Both the pattern transmission and distance to the focal plane affects the SNR of the information available to track. Forthis reason, the SNR is provided for comparisons across these signal-degrading parameters. The value reported is the mean of SNRs at the frequencies used for tracking computed from simulated images containing the pattern only and noise only. Figures 16A-D show the mean displacement error and standard deviation measured for in focus patterns of varying transmissions with a displacement of dx = [10, 0] . For WPC, the tracking precision begins to degrade for transmissions greater than α = 0.8 as seen in Figure 16A₁ or SNRs below 20 dB as seen in Figure 16B. However, because the mean measurements remain close to the ground truth values, the tracking accuracy can presumably be increased with repeated measures at the same pattern position. The pattern transmission has no effect on the MSC method in this scenario as seen in Figure 16C; and the error for this method is once again highly structured as seen in Figure 16D. At large SNRs, the error in the MSC method is significantly greater than that of the WPC method.

Focus

Figure 17 represents simulated square pattern with transmission α= 0.5 at different depths. From left to right in Figure 17, dz= {0, 5, 10, 15, and 20} μm. The test scenarios presented thus far have demonstrated the accuracy of the WPC-based tracking method for an in focus pattern. In structured illumination microscopy, however, the pattern can lie on a different plane from the specimen, and therefore is out of focus. The microscope PSF acts as a non-ideal low pass filter. As the pattern moves further out of focus, higher frequency components are attenuated to a greater degree, the pattern is smoothed, and the SNR decreases. Because higher frequency components are attenuated more than lower frequency components, the signal from high frequency components used in tracking will fall off faster than the signal from low frequencies. This may limit the usefulness of higher harmonics in the square wave pattern. Figure 17 shows the square pattern at different axial distances from the focal plane of the microscope.

Figures 18A-D show the effect of defocus on translation estimation for images, for example simulated with a 4OX, 0.65NA objective lens. Figure 18A illustrates translation error for WPC as a function of z position while Figure 18B illustrates translation error for WPC as a function of SNR. Similarly, Figure 18C illustrates mean translation error for MSC as a function of z position and Figure 18D as a function of SNR. The displacement of the square grid pattern from the focal plane ranges within dz = [0...30] μm (the analytical PSF is symmetric about dz = 0). The pattern is displaced laterally by dx = [10, 0]^τ . The data shows that the accuracy and precision of the displacement estimates fall off as distance from the focal plane increases. Tracking with WPC using four harmonics begins to deteriorate significantly above dz = 20 μm as seen in Figure 18A. Note, however, that monotonically increasing z position does not correspond to monotonically decreasing SNR. For example, the error for x translation at 14 μm is based on a SNR of about 2 dB. The results herein demonstrate the value of reporting the mean SNR at the frequencies used for pattern tracking. The WPC translation estimate degrades below 5 dB as seen in Figure 18B. The translation estimates of Figures 18C and 18D from MSC show fewer large errors, with consistent performance down to 0 dB SNR. At high SNR, MSC again under performs compared to WPC.

As a final comparison, Figures 19A and 19B show the effect of SNR on frequency-based orientation estimation. Figure 19A illustrates orientation error versus SNR for varying transmission, and Figure 19B illustrates error versus SNR for varying z position. In all cases, the pattern rotation is 0°. Though the focus test case covers a greater range of SNRs, the general trend is the same in both situations. The mean orientation estimate stays close to the ground truth value-within 0.015 pixel. This suggests that orientation estimates will be improved with a large number of samples, as is the case when tracking over a large number of video frames. The fact that the pattern orientation stays constant over the course of an experiment helps mitigate the effect of errors in individual orientation estimates.

Evaluation Using Real Images The test cases presented above provide upper bounds on tracking accuracy under the condition that no specimen is present. Real microscopy images were collected for example, to demonstrate the technique in practical application. The greatest limiting factor shown from experimental evaluation was manufacturing the semitransparent patterns. This section outlines a set of experiments that used real microscopy data from a proof-of-concept system.

Results obtained with the real system can be compared to a simulation that closely matches the experiment. It can be argued that deficiencies seen in the real system can be overcome by improving the pattern manufacturing process.

The real pattern was created using 1000 mesh Gilder transmission electron microscope (TEM) grids from Ted PeIIa, Inc. [16]. As mentioned previously, these grids are composed of repeating 25 μm square regions with

19 μm holes and 6 μm bars. To create a pattern, the grids are taped to a specimen holder or covering member, for example #0 (approximately 110 μm thick) glass slide cover and metal is deposited onto the glass, for example by using a thermal evaporator. The grids can then be removed from the specimen holder or covering member, leaving a metal layer on the specimen holder or covering member that is a negative image of the grid. Varying the thickness of the deposited metal can control optical transmission of the pattern. The metal deposition consists of 2 nm Cr and varying thicknesses of Au ranging between 10 and 30 nm.

This method of pattern manufacture yields highly variable results. Because the grids are fixed to the specimen holder or covering member along the outside edges, some grids warp during metal deposition, leaving unclear images of the grid. This leads to significant variation in optical transmission, even for metal deposits of the same thickness. These shortcomings, however, are purely a consequence of this particular procedure. Photolithography using precise masks could greatly improve the quality and repeatability of the pattern manufacture process, and would lead to obvious improvements.

Measuring Optical Transmission

Patterns were selected by visually assessing the clarity of the pattern images and by measuring optical transmission (the percentage of light transmitted through the metal layer). Knowing optical transmission of a pattern can help predict what sort of tracking accuracy should be obtainable. Pattern segments an image plane into two regions: a bright region, Rb where all light is transmitted, and a dark background region, Rd, where some constant fraction of iight is transmitted. For a constant transmission, α , the intensity levels measured in bright and dark regions can be represented in one embodiment as μRd = α μRb.

Transmission constant α can therefore be measured by segmenting a single image of the pattern into bright and dark regions, and by finding a ratio of light intensity measured in each region. For example, an expected value of a pixel measured by a CCD is given as μ (x, y) = k(x, y)E(x, y)A + μt(x; y)A, where k is a per-pixel gain from fixed-pattern and flat-field non-uniformities, E is the photoelectron count measured at a pixel, A is the ADC gain, and μt is an expected dark current value. After flat-field and dark-current correction, the set of pixels in one region provides an estimate of the intensity level in the region, ERA, that would be measured by an ideal CCD to be,

where | R | is the number of elements in a region R, and I _c is a flat-field and dark-current corrected image.

Using a single image to measure transmission, however, assumes that the image sensor has a linear response to intensity. Typically, CCDs usually have a very linear response, except at low intensities [15]. To account for the nonlinearity at low intensity levels, the equation relating intensity levels should be modified with a constant offset, as seen by: μRd = α μRb + β .

This modification assumes the camera has an affine response within the middle intensity ranges. A series of images of a pattern taken at different illumination intensities can verify that a camera has an affine response and can also estimate the optical transmission of the mask. This is completed by fitting the mean intensity values from the bight and dark regions in each image using the immediate equation above.

Pattern image segmentation is required to complete the transmission measurement. Segmentation divides the image plane into disjoint bright and dark regions. Ideally, there should be no overlap between regions, but there may be pixels that do not belong to either region. The intensity values of the pixels in each region have a normal distribution, thus K-means clustering provides the basis for an effective segmentation method [17]. K-means clustering, with k = 2, applied to pixel intensities in an image of the pattern divides the pixels into two clusters around two centroids with intensity values {IA, IB}- Each pixel may belong to a cluster whose centroid value is closest to its intensity value. The clustering is optimal in the sense that the sum of squared intensity differences between the centroid value and all cluster pixel intensity values is minimal over the two clusters. The initial clustering enables computing the expected values μRb and μRd and the standard deviation of intensities in the clusters, σ Rb and σ Rd .

Given these, values, a segmentation threshold for each region is provided by

{(x, y) e Rb I lc(x, y) > μRb - pσRb}

{(x, y) e Rd I lc(x, y) < μRd + pσRd}; where p is an empirically-determined constant that indicates how broad a region of intensities to include in each segmentation. Each segmented region may become subjected to a morphological close operation, which grows and then contracts region boundaries-to eliminate isolated pixels.

Figures 20A-C show a sample of images of a pattern taken with different lamp intensities. A sequence of pattern images obtained at different illumination levels is shown in Figure 2OA. Figure 2OB shows segmentation of the image plane with pixels R_d in medium gray and R_b in bright gray. Pixels belonging to neither region are in black. Figure 2OC shows mean dark intensity versus mean bright intensity obtained from 200 images at different intensity levels. The dotted line in Figure 2OC shows least squares regression. Figure 2OB shows the segmentation of the image plane based on the first of these images. Linear least-squares regression conforms that the camera sensor has an affine response over the middle of the intensity range, with a correlation coefficient, r2 = 0:998. This particular pattern has a transmission coefficient α = 0:580, which is among the lowest transmission coefficients obtained using the TEM grid manufacturing process. The offset β = -0:937 indicating image sensor is nearly linear. Matching Experiment to Simulation

Two experiments are presented here that evaluate two tracking methods using microscopy images of the square grid patterns. The first experiment evaluates tracking on the pattern grid absent a specimen. In one embodiment, capturing a series of images of a pattern moving with known displacement, on a prepared specimen holder or covering member is shown. The imaging system consists, for example, of an inverted Nikon Eclipse TE2000-E microscope with a 4OX 0.7NA objective lens and a Pulnix TM-671 OCL camera, which captures monochrome, 8-bit, 648 x 484 pixel² images at 120 fps. Microscope is focused on pattern which has a measured optical transmission of α = 0.46. The stage is moved in 10 μm increments along one lateral axis, with the Pulnix camera acquiring 5 images at each location. Each image is flat-field and dark-current corrected after capture. The MCL stage has a 100 μm range of motion along each axis (x, y, and z), but the total travel for this experiment is restricted to 80 μm to avoid intermittent misalignment issues reported at the far ends of the stage range. Phase-based and spatial correlation tracking is performed on 440x440 pixel² images to avoid introducing any bias from having more samples along one image axis. Figure 21 shows a selection of the images acquired during this experiment. The images are of a pattern being translated by a MCL nano-positioning stage. Figure 22 shows images of a pattern cropped from a single high-resolution image used to simulate stage motion.

A second set of images are used to simulate stage motion to closely match the experimental images. A high resolution image of the pattern is formed for example, using a Leica DM5000 B microscope with a 4OX, 0.85NA objective lens and a 1.2X tube lens for a 48X total magnification. Images are acquired for example, with a SPOT Flex FX1520 camera, configured to capture monochrome, 12-bit, 1024 x 1024 pixel² images. The high resolution image was scaled to match the intensity range and transmission of the experimentally- acquired images. Cropped regions of this image were selected to simulate translating the stage in 10 μm increments over 80 μm. Each image was corrupted with Poisson-distributed noise with ADC gain of 100 photoelectrons per count as measured, to create 5 simulated images at each location. Figure 24B shows a selection of simulated images. Table 1.5 illustrates RMS and standard error for tracking a pattern using WPC and MSC. Images are of a 25 μm square grid pattern in focus.

Table 1.5 above records the RMS tracking errors obtained on the experimental and simulated images using the WPC and MSC methods. Note that the bars in the pattern of Figure 21 are much narrower than the bars in simulated pattern images, for example Figures 10. This can be because metal diffuses around the TEM grid bars during the deposition process, narrowing the region that should remain completely transparent. The amount of diffusion is different for every pattern manufactured with this process, and the disparity will have a larger effect on the MSC method, as the grid in the modeled images will completely cover the grid in the observed images for a wider range of model parameters.

A second experiment investigated pattern tracking in the presence of a specimen. The pattern-prepared glass specimen holder or covering member is mounted on tissue sections prepared by biologists at the University of North Carolina at Chapel Hill (UNC). The specimen holder or covering member is placed pattern-side down, and secured with mounting media, which can leave a small gaps between the specimen and the pattern. Stained octopus muscle sections serve as specimens for this study. Figure 23 shows the octopus muscle and pattern at 10X magnification. In Figure 23, the metal deposit grid is mounted on a section of octopus muscle tissue. The pattern is clearly visible in the gaps near the center section. The edge of the pattern is near the bottom of the image. The pattern here is approximately 8μm above the tissue specimen, described below.

Image capture of experimental data follows a similar process to that used in the grid-only experiment above. A prepared slide is placed on a MCL nano-positioning stage and imaged with for example, an inverted Zeiss Axiocam microscope with a 4OX 0.7NA objective lens and a PulnixTM-6710CL camera. Microscope is focused on the specimen tissue, and pattern is approximately 8 μm above the specimen. The measured optical transmission of the pattern is α = 0:837. The stage moved in 10 μm increments along the x axis a total of 80 μm, with camera acquiring 5 images at each location. Each image was flat-field and dark-current corrected after acquisition. Figure 24A shows a selection of the images acquired during this experiment. For example, Figure 24A shows the experimental images acquired using a real pattern, and Figure 24B shows simulated images created to match the experimentally acquired images.

In Figure 24B, one embodiment of a simulated image was formed for example using a Leica DM5000 B microscope with a 4OX, 0.85 NA objective lens and a 1.2X tube lens for a 48X total magnification. 12-bit, 1024x1024 pixel ² images were acquired for example, using a SPOT Flex FX1520 camera. Pattern image is simulated as described earlier with a 40X, 0.65 NA objective.

Pattern is given an optical transmission of α = 0:837 and displaced 8 μm from the focal plane. The high resolution image was scaled to match the intensity range of the experimentally-acquired images, and multiplied by the simulated pattern image. Cropped regions of the modulated image were selected to simulate translating the stage in 10 μm increments over 80 μm (measured by the simulated pattern). Each image was corrupted with Poisson-distributed noise with ADC gain of 100 photoelectrons per count to create 5 simulated images at each location. Figure 24B shows a selection of simulated images.

Both sets of images can be tracked using, for example, MSC and WPC. All image frames can be registered to the first frame in the image sequence. Pattern spacing parameters in experimental images can be optimized for a fundamental frequency and four harmonics from images of the pattern alone. These four frequencies are used to perform the phase correlation. In table 1.6, RMS and standard error for tracking a pattern and specimen using WPC, MSC, and image registration are shown. Tracked images are of octopus tissue and a 25 μm square grid pattern, approximately 8 μm above the tissue.

RMS tracking error and standard error for the different tracking scenarios are enumerated in Table 1.6 above. Both methods perform worse on the octopus data set than the grid-only data set for example, with the error in WPC doubling and the error in MSC increasing about six fold. This data provides very important insight to the error characteristics obtained in simulation which match well with experimental data. Considering the x and y axis errors together, both WPC and MSC demonstrate approximately a 1.6 μm (7.1 pixel) RMS error in both the experimental data and the simulation data. The data indicates that simulations including corruption from a specimen provide a reasonable substitute for experimental data.

Creating Microscopy Mosaics

To close the chapter on Structured Illumination Microscopy Tracking, an example of WPC-based pattern tracking is given. In one embodiment, pattern is introduced to optical train of bright-field microscope to provide stage position information for non-motorized stages. Stage position information can be used to construct mosaic images from multiple frames of a microscopy video. A complete tracking application can enable microscopists to examine a specimen while a camera continuously acquires images. Stage position can be tracked, and the acquired images can be assembled and presented to the microscopist with proper spatial context. Assembled images can enable the examination of regions larger than a single FOV at high magnification. This provides a natural view of specimen structure, and may better enable microscopists to understand the relationships between different regions in a specimen.

When moving a stage micrometer by hand, it is difficult to restrict motion to less than half the pattern wavelength in each direction (for the TEM grid pattern, 12.5 μm) for every frame of an image sequence. In this case, WPC tracking is subject to an ambiguity of whole pattern distances. Image correlation can easily resolve this ambiguity for fixed specimens-image registration restricted to several whole-pattern displacements along each pattern axis provides a disambiguation of image pair alignments. WPC or MSC on the full- resolution image provides fine alignment.

Figure 25 shows an image mosaic assembled from 20 microscopy images using this approach. Images were formed, for example, on a Leica DM5000-B microscope with a 4OX, 0.85NA objective lens and a 1.2X tube lens. A SPOT Flex FX1520 camera acquires monochrome, 14-bit, 1024x1024 pixel images at 4 fps with 5 μs exposures. The stage was moved continuously by hand during image acquisition. During tracking, normalized cross-correlation on images reduced to 256 x 256 pixel displaced by [-3...3] pattern cycles disambiguates among possible translations reported by pattern tracking. The microscopy mosaic image was constructed by merging all images into a common reference frame, and then overwriting data with each new frame (no blending was performed). The resulting mosaic appears visually coherent, with some image seams visible upon close inspection. Improving pattern manufacturing processes could improve the tracking accuracy, as discussed next. This example demonstrates that the WPC tracking technique provides reasonably reliable information in some situations. One would not usually use structured illumination to create mosaics of the type seen in Figure 25. In fact, image registration on these images obtains a better alignment with fewer visible seams, and image blending techniques would improve the visual coherence of the mosaics even more [18]. However, WPC can additionally handle sparse specimens, moving specimens, and z position tracking.

In sum, a structured illumination approach for stage tracking in bright- field microscopy has been introduced, using theory and examples. A semitransparent pattern layer can be introduced to the microscope's light path, fixed to the specimen's frame of reference. The pattern transforms the stage tracking problem into a multiple layer image analysis problem. A search for magnitude peaks at pattern frequencies in the frequency domain can yield an estimate of pattern orientation that does not rely on knowing the pattern position. Two methods have been proposed to track the pattern. First, model- based spatial correlation (MSC) can determine pattern translation using an image formation model to simulate multiple images of the pattern, and optimizing model parameters to find the best match to an observed image. Second, weighted phase correlation (WPC) can compute a phase estimate for individual target frequencies in an image by first reweighting the image according to the target frequency's wavelength, and then computing a single component of the discrete Fourier transform (DFT). Once pattern's position has been estimated, a model-based rotation optimization can provide an improved estimate of the pattern orientation.

Accuracy of the orientation and tracking methods were examined in detail. For example, frequency-based orientation method provides estimates with a 0:25° RMS error for an in focus square grid pattern. Model-based orientation method improves this estimate to 0.066° RMS error once the position of the pattern is determined.

A different story can arise for pattern tracking, for example. WPC technique demonstrates very accurate tracking for axis-aligned square grid patterns without a specimen present. Lateral tracking accuracy remains around 0.010 pixel (2.3 nm for the simulated microscope) RMS error in this case, falling to 0.057 pixel for arbitrary orientations. This accuracy is maintained down to approximately 18 dB SNR due to pattern transmission and below 10 dB due to defocus effects. The WPC method relies on a phase estimate from individual pattern frequencies. Point sampling in the frequency domain requires infinite support in the spatial domain, and WPC attempts to obtain this by redistributing and reweighting the signal to analyze the target frequency. In the presence of uncorrelated noise, the WPC method maintains subpixel tracking accuracy. In the presence of correlated noise, however, the phase estimate may become unreliable, and tracking accuracy diminishes. Nevertheless, there are some practical imaging scenarios in which WPC may excel. For example, multimode microscopes capture fluorescence and bright-field images in quick succession [19]. For observations of fluorescent specimens that do not show up well in bright-field, one could use the bright-field channel to obtain tracking information from structured illumination. Because fluorescence microscopy uses epi- illumination, placing the pattern below the specimen for example, would cause no interference with the fluorescent imaging. In comparison, MSC struggles to provide 0.25 pixel RMS error in the absence of rotation, though it seems to improve remarkably for many non-axis- aligned patterns. Uniform accuracy is maintained down to below 10 dB SNR due to pattern transmission and to 0 dB due to defocus effects.

The error obtained with both methods increases significantly on real microscopy images of pattern grid. Tracking can degrade even more in the presence of a specimen. However, errors for experimental and simulated images with a specimen present agree quite well. For tracking contiguous, fixed specimens, one could use an image registration technique if in the presence of sufficient specimen contrast, as image registration may fare better than either WPC or MSC. However, these methods can provide something that image registration cannot-they enable tracking with a fixed frame of reference in the presence of sparse specimen information and for moving specimens.

Structured illumination provides the additional ability to extend tracking into 3D.

Pattern selection is discussed further below.

V. Pattern design for Structured Illumination Microscopy In Section IV above, weighted phase correlation (WPC) and model- based spatial correlation (MSC) were developed as techniques that can provide lateral tracking of a semitransparent pattern layer in structured illumination microscopy. The pattern used in this evaluation was obtained from a TEM grid which produced dark square \holes" separated by thin bright "bars". This pattern was chosen because the grids were readily available to purchase from a microscopy supply company and it was relatively simple for a materials scientist to deposit a metal image onto glass cover slips. This provided a working system with which to demonstrate the concept of structured illumination microscopy. The TEM grid pattern, however, is not the optimal pattern for structured illumination microscopy, and pattern design for structured illumination will be discussed. The end goal in pattern selection includes gaining an understanding of how choices made in pattern design affect the ability to recover tracking information from the pattern layer. Specifically, how the choice of pattern influences the ability to track the stage independently from the specimen, determine focus position, and disambiguate among possible stage positions. Patterns may be formed by printing, etching, or photolithography and are not limited to such methods.

The following design goals are generally applicable to stage tracking in a broad range of microscopy experiments:

1. The pattern should enable tracking laterally and axially (in X, Y, and Z) independent of specimen motion.

2. The pattern should enable tracking over large frame-to-frame displacements.

3. The pattern should enable determining the absolute stage position.

4. The pattern should be easily removable from the acquired images through image processing.

Note that there is an important distinction between goals 2 and 3. The former goal refers to tracking the stage relative to some fixed reference position. The latter goal refers to determining absolute stage position from a single image-a slide-based absolute positioning system. Pattern design necessarily involves some trade-offs. For example, both

WPC and MSC can depend on finding the displacement of a regular pattern layer present in an image. The resolution of these methods increases as that pattern matching is based on higher frequency components, as seen in the equation below:

Epx = Eφ

36Of

But, to enable tracking over large frame-to-frame displacements the pattern should have a large wavelength, favoring low frequencies. Another consideration that favors low frequencies is that as the pattern moves farther out of focus, the high frequency components are attenuated more and therefore reach the noise floor first.

In order to track the stage independent of the specimen the pattern signal should dominate the specimen signal. This requires high-contrast, low transparency patterns. Patterns that are more transparent, however, enable removing the pattern from the specimen images with less information loss.

To enable axial (Z direction or focus) tracking, the pattern should contain multiple frequency components. The more frequency components present, the more terms can be used to determine focus. Choosing which frequencies to include, however, requires considering the frequencies present in the specimen. If the specimen has high magnitude at many frequencies, it may be difficult to find an appropriate set of frequencies to use as a pattern.

Additionally, the microscope point-spread function (PSF) may attenuate a large proportion of high frequencies too much to be useful.

Pattern Design

While it may be ideal to be able to design a single pattern that could provide optimal tracking in all scenarios, the reality is that tracking performance is largely dependent on the imaging application. For example, the size of pattern in images depends on the magnification provided by microscope and the size of camera's image sensor. Furthermore, some specimens have strong magnitudes at particular frequencies which interfere with the signal provided by the pattern. Coordinate System

One can specify wavelengths and frequencies in object or image coordinates, depending on what parameters are being considered. For example, in one embodiment, the relationship between a wavelength in object coordinates Lo, and the corresponding wavelength in image coordinates

(measured in pixels), Li, is

S

A spatial frequency in object coordinates (measured in cycles per μm, for example), f ₀ , can be represented as:

and a spatial frequency in image coordinates (measured in cycles per pixel), fi, can be represented as: f - ! ₌ ^s ₌ ^sfo 1 L₁ ML₀ M ^'

Pattern Spacing

Imaging application influences the selection of pattern's physical size. For example, in one embodiment the lower limit in pattern size is determined by three factors. First, the numerical aperture (NA) of the objective lens, second the spacing of the sensor elements in the CCD used to record images, and third is shot noise. An Abbe limit for example, relates objective NA to the smallest spacing in a resolvable pattern represented by: λ d =-

2NA

The Wittaker-Shannon sampling theorem contained within the References section states that the Nyquist frequency, f _N , is the maximum frequency represented by a sampling rate f _s [11] as:

f - ^fs In other words, accurately reconstructing a signal containing a specific maximum frequency requires sampling the signal at twice that frequency.

The size of the smallest pattern wavelength that can be resolved in an image taken with a CCD can be determined by twice the distance between

2s sensor elements, — , where M is the magnification of the objective lens and s

M represents the size of a CCD sensor element.

The noise in acquired images can be reduced using a Gaussian filter with standard deviation equal to the standard deviation of the noise expected in the image, σ n [13]. Applying a low-pass filter to structured illumination images can diminish the signal from patterns, especially those with short wavelengths. In the spatial domain, a Gaussian filter kernel with standard deviation σ can have 98% of its energy contained in a width of 5σ (or radius 2.5 σ ). So for example, an adequate guideline may be to choose pattern wavelengths that exceed at least twice the noise level to avoid excessive smoothing of the pattern signal during filtering.

Taken together, the Abbe limit, Nyquist rate, and noise level provide a lower bound for an observable pattern wavelength:

The upper limit to pattern spacing can depend on the field of view (FOV). Higher magnification applications require smaller pattern structures to include more than one pattern instance within a single FOV. For example, in patterns with wavelengths greater than one FOV, the WPC method may not have enough information to determine how to weight each component in its correlation computation. Though the mathematics still holds, the accuracy of the phase estimates is severely diminished. Additionally, the structure of the specimen places constraints on pattern spacing. Takita et al. note that natural images have a predominance of low frequency information [10], so there may be significant interference between pattern and specimen signals in this region. Although not as firm as the lower bound, the upper bound for a usable pattern size is on the order of the field of view along one axis of a microscope image:

where N is the number of pixels along an image axis.

Choosing Optimal Pattern Frequencies Assuming that pattern and specimen move independently, pattern can provide the only information in an image for tracking with structured illumination microscopy. To provide the most accurate tracking, maximize the information present in the pattern (perhaps the signal) with respect to the other factors in the image (for example, the noise). In this sense, noise can include both camera noise and specimen, as both of these elements interfere with pattern tracking.

The image that a microscope forms of a planar object is the convolution of the image predicted by geometrical optics with a slice of the three- dimensional (3D) point spread function (PSF) determined by where the object is placed in front of the objective lens. Equivalently, in the frequency domain, the magnitudes of the frequencies in the object are attenuated by the Fourier transform of the microscope's PSF by:

I_i(fe,fy) = h(fe,fy)I_g(&,fy) , where l(fx, fy ) is the Fourier transform of l(x, y) and h is the Fourier transform of the normalized PSF, also known as the optical transfer function (OTF).

Analysis of microscope OTF enables choosing of optimal pattern frequencies for a tracking application. The OTF is the frequency domain representation of the normalized point-spread function it determines by how much pattern frequencies are attenuated. So, maximizing the SNR for a pattern involves choosing the pattern that, when multiplied by the OTF, has a maximum magnitude at the pattern frequencies. The selection of possible patterns and OTFs from all possibilities is driven by the constraints of the experimental parameters. The critical experimental parameters can include for example, objective lens characteristics, image sensor size, the expected axial distance between specimen and pattern, and the composition of the specimen itself. Figures 26A and 26B show a frequency-focus slice of magnitudes in the OTF of a simulated 4OX, 0.65NA objective lens. In Figure 26A, one slice of the OTF displayed as an image is shown. Magnitudes are mapped with a logarithmic scale to make all magnitudes visible. Bright values indicate high magnitudes. Figure 26B discussed further is the OTF as a height map with scaling.

Here focal depths range within z = [0...25] μm and frequencies range within f = [0...1/2] cycles per pixel. For this lens and sensor pair, the Abbe frequency limit is above the Nyquist frequency. This figure illustrates why one cannot select a set of pattern frequencies that provide accurate tracking capabilities when the pattern is in focus, and expect tracking to perform uniformly throughout the focus range. Different frequencies can be attenuated by different amounts at each focal depth, and this attenuation is not monotonic across frequency or depth. Let P(x) be a pattern defined by a set of frequencies f := {f ; | i = 0,

1 ,...n}. Let P(f) be the Fourier transform of P(x). The magnitude at frequency f _j can be given by|P(f_i) . Let H(x; z) be the normalized PSF of an objective lens at a focal plane z. Here, the semicolon specifies that the focus axis is treated differently from the lateral axis on which the pattern resides. Specifically, H(f; z) is the Fourier transform of a slice of the PSF, one slice of the OTF; the magnitude at frequency f _j is given by |H(f_j;z)| .

The total magnitude at all frequencies in a pattern at a distance, z, from the focal plane is given by

The equation above provides a scoring function for how strongly a particular pattern with frequencies f_j is represented at a particular focal depth, and this indicates how well this pattern will trackable compared to other patterns with different frequencies. There is usually a known range of distances from the focal plane to the pattern in a given experiment. This distance in one embodiment, could be fixed, or it could stay centered around a value with some variation, or it could span a whole range with equal probability. For this reason it is sensible to introduce a weighting function W(z) that weights the importance of focal distances in an experiment. This can be thought of as a probability density function for axial stage positions.

Combining the pattern scoring function for individual focal depths and the weighting function for the importance of focal depths, a metric for how well a particular pattern is suited to tracking in an experiment is given by:

S(P) lPCf₁) IHCf₁Jz .

Note that the weighting function in the equation above is multiplied by all frequencies in the pattern, and the result is summed across all focal depths. Rearranging terms therefore provides a function that specifies the score for all frequencies in the OTF, irrespective of the magnitude in a particular pattern:

S(f) =∑W(z)|H(f;z)| .

Z The frequency selection thus far has not taken into account noise present in the system. Camera noise can have a constant magnitude at all frequencies. In the equation above, all frequencies are treated equally. Any frequency that is heavily attenuated by the OTF shall receive a low score, and so provides the set of best frequencies to choose. Spectral analysis of the specimen to be observed provides a way to consider the tracking interference caused by its presence in images. Because the specimen can have any arbitrary orientation, the magnitudes at frequencies at all orientations should be considered. Given a specimen image, l(x, y), and its Fourier transform in polar coordinates, l(f; θ ), the radially-summed magnitude provides a measure of specimen interference:

Q(f) =∑|l(f,θ)| . θ

To adjust this to the same scale as the axial weighting function, it is sensible to express the magnitudes on a logarithmic scale, and shift the scores to fall within Q(f) G [0...1]. This weighting term can now apply an additional constraint to the inverse of the specimen magnitude at all frequencies, Q'(f) = max(Q(f)) - Q(f), such that: S(f) = Q'(f)∑W(z)|H(f;z) .

Z

The frequency weighting function can also be truncated to enforce the minimum and maximum pattern frequencies.

Now that a frequency scoring criterion has been established, the final task is to use it to select which frequencies to include in the pattern. Choosing the n frequencies with maximum score maximizes the equation immediately above, but these frequencies are likely to be bunched around whichever frequency had the maximum score. For this reason it makes sense to apply some additional constraints to the frequency selection. For example, selection of frequencies that have local maxima in the scoring function distributes the pattern frequencies throughout frequency space, but uses the best one in each area. Or, selection of the frequency with maximum score and some of its harmonics ensures the pattern is periodic over a small window.

In summary, to determine the optimal frequencies to use in a pattern: 1. Select an objective lens and compute its OTF.

2. Construct a focus weighting function, W(z) that specifies the importance of different axial distances from the pattern.

3. Construct a frequency weighting function, Q(f) that specifies the importance of different frequencies, using the parameters of the imaging system and the Fourier transform of an image of the specimen.

4. Multiply the OTF by the focus and specimen weighting functions.

5. Sum the weighted OTF along the focus axis to determine the frequencies that will provide the best tracking signal.

6. Select the n frequencies that maximize this signal. Figures 27A-D provide some concrete examples of weighting functions applied to the OTF of Figure 26. In Figure 27A, the only constraints applied are the minimum and maximum frequency thresholds. The high frequency cutoff occurs because of the noise filtering requirement, where σ n = 1 :6, as determined for example by a Pulnix camera. For this scenario, the predominant feature of the frequency scoring function is due to the attenuation of higher frequencies by the OTF. Low frequency patterns will provide the best tracking signal. The vertical dashed lines in the frequency score plot indicate the four local maxima with highest frequency, but in this case there is little to distinguish these points from any others.

In Figure 27B, the focus range has been restricted to fall within dz =

[10...15] μm of the pattern. The frequency score again shows a preference for low frequency patterns, with local maxima at 0.025 and 0.05 cycles per sample that would be good options for patterns that have frequency spread throughout frequency space.

In Figure 27c, the focus range is restricted within dz = [0...15] μm and the frequency weighting has been constrained by the inverse spectrum of a frog brain tissue specimen. Because the tissue has a predominance of low frequencies, the optimal frequencies to include in the pattern are pushed out to around 0:2 cycles per sample.

In Figure 27D, the frequency weighting is the same as above, but the focus weighting has been provided by a normal distribution of stage positions, with dz ~ N(17 μm, 2 μm). In this case, a mixture of low frequencies and high frequencies could provide the best tracking pattern.

The Tracking Noise Floor

Once the pattern frequencies have been chosen, the question remains of how well pattern will be trackable in an experiment. An answer to this question is provided by analyzing the signal-to-noise ratio (SNR) of the system, where once again the signal is provided by the pattern and the noise is provided by camera noise and the specimen.

The SNR is the ratio of information in the signal to interference in the noise. The information present in a signal is given by the magnitude at the frequencies in the signal. The relationship between signal amplitude in the spatial domain and its magnitude in the frequency domain is:

where N is the number of samples in the signal. The maximum amplitude in a signal is determined by the dynamic range required to represent it. For digital images, dynamic range may be measured by bit depth, where the dynamic range of an image, l(x, y), is

"Oτ*~" _{b = ]og} fmaxCIOyO)^" ^ min(I(x,y))

At maximum amplitude, a signal contains both the highest and lowest values allowed by the dynamic range. The maximum amplitude for a signal represented by b bits is then represented as: A ₌i_ _ ob-l

and the maximum magnitude at any frequency is

M_m Iϊl_aα_xX 2^b~2N .

Given a pattern with dynamic range b bits, a specimen with dynamic range b _s bits, and noise of standard deviation σ n, the signal to noise ratio for saturated pattern and specimen is:

SNRff z) - ^{2V2 NH}(^{f> z})

2^{bs "2}N + M_n where Mn is the mean magnitude of the noise, provided by Mn = σ nVN . his is a worst-case estimate of SNR, as it assumes that the specimen has maximum dynamic range at all frequencies. It is preferable instead to compute the SNR only at the frequencies in the pattern.

When the SNR falls below 1 (or 0 dB), the pattern signal can be said to be below the noise floor. In practice, the signal from a pattern is required to remain some amount above the noise floor.

Increasing the Signal-to-Noise Ratio

In one embodiment, increasing the SNR requires either increasing the pattern signal or decreasing the specimen signal. The dynamic range required to record the bright and dark regions in a specimen is fixed-that is, a specimen that covers 4 bits of dynamic range to record detail in its brightest and darkest regions on for example, an 8-bit camera will also cover 4 bits of dynamic range on a 12-bit camera. For a specific specimen, an increase in pattern SNR can therefore be obtained only by increasing the pattern's dynamic range. The pattern's optical transmission can become limited to the range α = [0...1], and every halving of the minimum optical transmission increases the pattern dynamic range by one bit and consequently decreases the dynamic range remaining for representing the specimen by one bit. Increasing the camera's dynamic range and decreasing the optical transmission of the pattern can provide a practical method to increase tracking SNR without affecting the dynamic range used to record specimen information.

Another factor that can affect pattern's SNR is the distance of the pattern from the focal plane. At larger distances from the focal plane, more light diffracts from the bright regions of the pattern into the dark regions, lowering the pattern's effective dynamic range. Note that the optical transmission of the pattern remains the same, but the ratio of light in the brightest and darkest regions of the pattern diminishes. Therefore another practical method to increase the pattern SNR is to decrease the distance between the specimen and pattern.

Determining Focus (Axial, Z Tracking)

A pattern used in structured illumination microscopy can provide axial as well as lateral tracking. The method relies on analyzing the defocus of the pattern in observed images. In a focus model the PSF is symmetric about the optimal focus plane (z = 0) and is usually radially symmetric about the optical axis. The microscope stage can moved axially (in z) to align a plane of interest in the specimen with the focal plane. As the stage moves, different parts of the specimen can become in focus at the image plane while other parts move out of focus. As discussed earlier, frequencies in the pattern can be attenuated by the OTF dependent on the distance from the focal plane.

Knowledge of the pattern and microscope PSF leads to a method for axial tracking of the microscope stage. In one embodiment, the focus model predicts the attenuation of pattern components expected at different focal distances. The model predictions can be compared to experimental images to obtain an estimate of the pattern distance from the focal plane-this estimate provides axial tracking for the microscope stage. Because observations are made with the specimen in focus, knowing on which side of the specimen the pattern is placed resolves the ambiguity from the axial symmetry of the PSF.

For example, let f_j = (r_j ; θ_j ) be a set of frequencies known to be present in a pattern, specified in polar coordinates where T₁ is the distance from the center of frequency space and O₁ is the orientation at which the component appears. The magnitudes observed in an image l(x, y) are given by |l(iϊ,θ_j with the transform to polar coordinates. The magnitudes predicted by a model pattern image P(x, y; dz) for a distance from the focal plane dz are given by P(r_j,θ_j;dz . Given an observed image and a model of the microscope and pattern, focus estimation seeks to minimize an objective function, for example by:

E(dz) = ∑(1(_ΑA) P(r_i5θ_i;dz)

where the distance between observed and modeled magnitudes is given by the

I norm. Focus estimate can be applied to any number of frequencies present in the pattern, but that number may be limited depending on pattern design. Due to the non-monotonic nature of the OTF, focus estimation objective function can be characterized by many local minima over practical focus ranges. Arriving at the global minimum and a reliable estimate for z position relies on either a search at a fine sampling of the objective function or restricting the search range based on experiment parameters.

The former option involves a practical balance between sampling resolution and performance, as evaluating the PSF requires significantly more processing for larger axial displacements. This performance hit can be mitigated by first conducting a coarse search through focus, applying optimization around all local minima, and selecting the best one. This is the approach adopted here, where local minima within 10% of the smallest value in the range of focus scores are considered.

The latter option is also practical in many imaging situations. One can obtain an initial estimate of the distance between the specimen and the pattern using a stage micrometer. Once an initial focus range is established, one can restrict the search space by assuming that focus will change little from frame to frame. That is, the focus search can be conducted within a bounded range around the prior focus estimate. With a suitably-constrained search range, optimization of focus objective function can proceed with a 1 D minimization technique, such as golden section search [7].

As an example, Figures 28A and 28B demonstrate focus estimation for a simulated sinusoidal pattern using four frequencies along each pattern axis. The pattern in Figure 28A has transmission α = 0:25 and is placed 6 μm from the focal plane. Images are simulated for example, using a 4OX, 0.65NA objective lens, 9x9 μm CCD image sensor, an ADC gain of 100 photoelectrons per pixel and Poisson-distributed noise. Figure 28B shows the objective function of the equation immediately above evaluated for displacements of dz = [0...25] μm with the observed image in Figure 28A. With this sampling for example, the minimum occurs at 6 μm, but another local minimum occurs at 16 μm. So, it appears that accurate z tracking depends not only on there being sufficient pattern SNR, but also upon ability to discriminate among different regions of the focus estimation function. In this case for example, optimization in the regions near the local minima resolves a focus estimate of 6.01 μm. Focus measure can depend on knowing the orientation of the pattern in the image. A method for estimating pattern orientation was discussed earlier. The focus measure does not depend on the lateral position of the pattern, but rather from the Fourier shift theorem, a translation in the spatial domain affects only the phase in frequency domain. The WPC tracking method does not depend on focus, but MSC does. These considerations suggest the natural order for practically obtaining stage tracking information. First, an estimate of orientation of the pattern in the frequency domain from the maximum magnitudes at each pattern frequency is made. Next, an estimate of axial (Z) position of pattern is made by comparing the magnitudes of the pattern frequencies to magnitudes predicted by the focus model. An Estimate the lateral (X₁Y) position of the pattern can be made using either WPC or MSC. Finally, the orientation estimate can be refined using MSC, and these steps may be repeated until convergence. Evaluation Using Simulated Images

Evaluating a test scenario begins with choosing appropriate system parameters, for example, the pattern design, objective lens, camera parameters, and specimen. Test images can then simulate for a variety of stage translations with noise added. Stage position is estimated for each translation by averaging the translation obtained by the tracking algorithm for a number of different noisy images. Given the errors for each translation the root mean squared (RMS) error is computed for the test scenario.

Pattern Spacing

To evaluate the effect of pattern spacing on tracking accuracy, tracking is performed on simulated images of single frequency sinusoidal patterns. The

10 512 wavelength of the patterns ranges within L_n = [ ... ]. For each

256 256 wavelength, the pattern can be translated along the x axis between dx =

[ — ... — — ] at lO even spacings. For each image, the pattern displacement is

measured with WPC and MSC in reference to an untranslated pattern image. The pattern is given a transmission of α = 0:25, placed dz = 5 μm from the focal plane, and imaged with for example, a 40X, 0:65NA objective lens to form

8 bit, 256x256 pixel images. Each image can be corrupted by Poisson- distributed noise modeled by 100 photoelectron/pixel ADC gain.

Figure 29A shows the RMS tracking error as a function of pattern spacing in x and y for each tracking method as a function of pattern wavelength. It is immediately apparent that WPC tracking deteriorates significantly for wavelengths longer than %th the FOV. Figure 29B shows the same data as Figure 29A, but zoomed in to the region marked by the dashed lines. MSC maintains consistent tracking accuracy for wavelengths up to 1 :2 times the field of view. Figure 29B shows a subset of the data in Figure 29A.

This view shows that the RMS error for MSC remains below 0:02 pixels in most cases. Tracking Without Specimen

Structured illumination patterns containing low frequencies are optimal in the absence of a specimen. The tracking test cases from Section IV above were rerun for an optimized pattern consisting of orthogonal sinusoidal waves.

Table 2.0contains a comparison of the RMS error for tracking square grid and sinusoidal patterns. Bold values indicated the best accuracy obtained for a particular test case considering both x and y accuracy.

Table 2.0compares the RMS errors obtained for tracking the grid pattern from Chapter 5 and optimized sinusoidal patterns using WPC and MSC. The test cases refer to the orientation and translation evaluations performed Section IV, and the data for tracking the grid pattern comes from this section, repeated here for ease of comparison. In this evaluation, all simulated patterns have transmissions of α = 0.5, are imaged in focus with for example, a 4OX, 0.65NA objective lens to obtain 8-bit, 256 x 256 pixel images. The fundamental plus three higher harmonic frequencies are used for tracking the square grid pattern; the sinusoidal pattern contains the four frequencies with maximum score along each axis. The bold values in this table represent the pattern and tracker combination that obtained the best accuracy for a particular test case, 2 evaluated by smallest RMS magnitude, E =y RMS_x + RMS . In all cases the optimized sinusoidal patterns tracked by MSC showed the least error, remaining below 0.005 pixel (1.1 nm).

The objective function for the sinusoidal patterns used in MSC is well- suited to gradient-based optimization, where the direction in which to search for a minimum is obtained from an estimate of the objective function's gradient. Gradient-based methods are known to converge faster than Nelder-Mead simplex-based optimization when applied to the right problem [7]. Applied to the square grid pattern, gradient-based optimization converges more slowly and to a less accurate solution compared to simplex-based search. Applied to the sinusoidal pattern, gradient-based optimization converges more quickly and to a more accurate solution than simplex-based search. Therefore, switching from the grid pattern to patterns composed of relatively few sinusoids enables using a faster optimization strategy that obtains more accurate results.

Tracking Moving Specimens

Section IV demonstrated that tracking a pattern layer using MSC and WPC was possible in the presence of a specimen that moved in concert with the pattern. This situation is common for example, in histology studies, but in many such cases specimen-based image registration techniques can also provide the tracking data. An obvious exception is when there is not enough specimen information to provide reliable image registration, for example in gaps or sparse sections within a specimen or between sections of tissue.

The full power of structured illumination microscopy includes the ability to track a pattern layer independent of the specimen motion. This expands stage tracking capability from observations of stationary specimens to experiments involving live, moving specimens.

Sparse Specimens Figure 3OA shows a simulated image of a specimen with sparsely- distributed contrast. Figure 3OB shows the pattern composed of three optimal frequencies and one low frequency. Figure 3OC shows the SNR estimated for tracking the pattern in 3OB in the presence of specimen in 3OA as a function of pattern distance from focal plane.

The image in 3OA represents what a microscopist may see when looking at a cluster of opaque beads, for example in a bead diffusion experiment, or at a colony of bacteria. Optimal pattern selection on this specimen image, as discussed earlier, reveals that the best frequencies to use to track a pattern in the region dz > 5 μm from the plane appear around 0.18 cycles per pixel. Because patterns at this frequency do not allow for much stage displacement between image frames, the low frequency with the highest score was also selected to be included in this pattern. Figure 3OB shows the pattern generated with these optimal frequencies. Figure 3OC shows the SNR for these frequencies as a function of pattern z depth, computed at the pattern frequencies in the presence of the specimen and Poisson-distributed noise with an ADC gain of 100 photoelectrons per count. The legend on this graph displays pattern frequencies expressed as cycles per image.

Tracking the pattern in the presence of a moving specimen was simulated by moving the pattern and keeping the specimen fixed. To maximize the pattern SNR further from the focal plane, the pattern can be given a minimum transmission of α = 0. The pattern was simulated at different focal depths while keeping the specimen in focus. At each focal depth, the pattern

L. was displaced along one axis within dx = max L max .where h_max is the

maximum wavelength in the pattern.

Figure 31 A shows RMS tracking errors for tracking the optimized pattern at different distances from the in focus specimen. The RMS tracking error is in the presence of a specimen with sparse contrast. For this data set, MSC obtains tracking errors below 0.03 pixel for z positions less than 5 μm, outperforming WPC in this range. The tracking error for MSC remains below 0.5 pixel for z positions less than 15 μm. WPC outperforms MSC above 5μm, maintaining an error below 0.5 pixel to almost 20 μm from the focal plane. This represents a distance of 8.5 times the depth of field for the simulated lens.

Figure 31 A shows RMS tracking errors for tracking the optimized pattern at different distances from the in focus specimen. From this plot, the tracking methods maintain an accuracy of better than 0.5 pixel down to the noise floor (0 dB) for WPC and down to 5 dB for MSC.

Frog Brain Tissue To simulate tracking in the presence of specimens with densely- distributed contrast, test images are create using single real microscopy images that remain stationary multiplied by simulated pattern images. The microscopy images come from the Burmeister lab at the University of North Carolina at Chapel Hill (UNC) Department of Biology. Scientists Mangiamele and Burmeister used radioactive markers to localize the expression of immediate early genes that signal neural activity in tungara frog (Physalaemus pustulosus) brains [20]. Figure 32A shows a section of frog brain tissue; the gray regions are neurons and the dark spots are silver grains from photographic emulsion used to localize the radioactive markers. Although the brain tissue image remains stationary during these tests, moving the pattern simulates the effect of the specimen and pattern layers moving independently. Figure 32B shows the pattern composed of four optimal frequencies including one low frequency.

Figure 32C is the SNR estimated for the pattern in Figure 32B in the presence of the specimen in 32A as a function of pattern distance from the focal plane. Figure 33A shows RMS tracking errors for tracking the optimized pattern in the presence of the tissue at different distances from the in focus specimen. For this data set, MSC obtains tracking errors below 0.5 pixel for z positions less than 4 μm, outperforming WPC. This provides subpixel tracking for patterns that remain within one depth of field of the specimen. The tracking error for both MSC and WPC becomes greater than 1 pixel for axial distances above 5 μm.

Figure 33B shows the same RMS tracking errors as a function of mean pattern SNR, considering the frequencies in Figure 32A. From this plot, MSC maintains an accuracy of better than 0.5 pixel (0.1123 μm) almost down to 10 dB. Note that this tracking accuracy is well below the Abbe resolution limit of 0.423 μm for this lens. The tracking deterioration for axial displacements above 5 μm is well predicted by Figure 32C. This example further demonstrates why the tracking in Section IV suffered for the square grid patterns in the presence of specimens. The square grid patterns are composed of a fundamental frequency and harmonics at decreasing magnitudes. At low frequencies the specimen provides the greatest interference, diminishing the tracking accuracy for frequencies in the pattern with the largest magnitude. At higher frequencies, the magnitudes in the pattern are too low to rise above the camera noise.

Estimating Focus Focus estimation using the approach described earlier was evaluated for simulated experiments both with and without a specimen. The axial tracking range is limited by the depth of field of the objective lens-as the pattern moves farther outside the depth of field for the lens, the signal diminishes to a level that cannot be tracked reliably. The depth of field for a microscope system is represented as: λn ne d = ^• =^■ + ^■

NA² MNA ^' where λ is the wavelength of light used in the observation, n is the index of refraction of the objective immersion medium, and e is the resolvable distance of the image sensor-three times the physical size of the sensor elements. The embodiment discussed here assumes illumination with 550 nm green light, the use of a dry objective (n = 1.0), and a CCD sensor element size of 9 x 9 μm.

A pattern with minimum transmission α = 0 can be simulated at different depths from the focal plane using for example, a 4OX, 0.65NA objective lens, and Poisson-distributed shot noise modeled with 100 photoelectron per count ADC gain. This optical system has a depth of field of 2.34 μm. Axial distances to the pattern were within the range dz = [0...30] μm and the initial search range is restricted to dz = [-3...3] μm, with the assumption that axial position changes slowly from frame to frame. For each focus position, five images are generated with different noise, and the mean estimated z position is recorded.

Figure 34A shows mean focus estimation errors for simulated images of a sinusoidal pattern with four different frequency components along each pattern axis. The mean z position error as a function of axial distance to the focal plane is given in Figure 34A. The pattern frequencies were selected using the same criteria used for lateral tracking. In the absence of a specimen, low frequencies are preferred. The focus estimation error remains below 0.5 μm to at least 12.5 times the depth of field with the exception of when the pattern is completely in focus. The RMS focus error is 0.266 μm with a standard error of 0.244 μm. Figure 34B shows the effect of mean pattern SNR on the axial tracking, demonstrating that tracking persists down to the noise floor. Figure 34B is the mean z position error as a function of mean pattern SNR measured at pattern frequencies.

Figures 35A and B show the z tracking error for a simulated sinusoidal pattern multiplied by an image of an in focus tungara frog brain section. Figure 35A is the mean focus error versus z position and 35B is versus SNR. The frequencies for this pattern were selected using the optimal lateral tracking criteria, with the additional constraint that a wide range of frequencies are chosen which amounted to selecting frequencies with local maximum scores that were not among the absolute maximum scores. The tracking search range is bounded within [-3...3] μm of the ground truth value, applying the assumption that z position changes slowly from frame to frame. Adequate tracking accuracy occurs when the pattern is within dz = [0...15] μm from the focal plane, maintaining less than 0.5 μm error in all cases except at 10 μm, where there is a temporary drop in the SNR. The range extends at least 6 times the depth of field for the lens. In this range, the mean pattern SNR spans [35...5] dB. At lower SNRs the tracking accuracy deteriorates significantly. In contrast, lateral (x and y) tracking for the sinusoidal pattern in the presence of the frog tissue has a RMS error of 0.113 μm for SNRs above 10 dB. The theoretical axial resolution of a bright-field microscope is half the lateral resolution, so some degree of diminished accuracy is to be expected [21]. The reduced lateral tracking accuracy corresponds well with the theoretical reduced axial resolution for the bright-field microscope for example.

As a comparison, Figures 36A and B show the z tracking error for a non- optimized sinusoidal pattern in the presence of the frog brain section. Figure 36A illustrates mean focus error versus z position and Figure 36B is versus SNR. The ranges in these graphs are selected to match those in Figures 35A and B. Note that z distances only up to 20 μm are evaluated in Figures 36 A and B. Frequencies for this pattern were chosen arbitrarily, and the impact on tracking error is apparent-tracking error is less than 1 μm in the range dz = [2...8 μm]. The pattern SNR drops below the noise floor when the pattern is greater than 10 μm from the focal plane.

Increasing the z tracking range using structured illumination microscopy requires increasing the depth of field of the objective lens. Lower magnification lenses usually have smaller NAs, and consequently a significantly larger depth of field, as there is a dual effect on the depth of field. Intuitively, at lower magnifications and NAs, the PSF projects onto a smaller region of the image sensor, so out of focus elements are not smoothed as much over the image plane.

Figures 37A and 37B show the z estimation error for simulation experiments conducted for example, with a 2OX, 0.40NA objective lens, which has a depth of field of 6.81 μm. Focus error is plotted against z position in Figure 37A and against SNR in Figure 37B. The axial distance to the pattern ranges within dz = [0...30] μm and the initial focus search is constrained to within [-3...3] μm of the true focus position. The objective's increased depth of field enables obtaining z tracking estimates over a larger range. The focus estimate is inaccurate when the pattern is placed close to the specimen because the pattern appearance changes little in this region, and fitting the focus objective function is therefore particularly prone to error. The tracking estimate becomes reliable when the pattern is about 1 μm from the specimen. The estimate has less than 0:5 μm error within the range dz = [1...30] μm, and the mean pattern SNR ranges from [35...3] dB over this range.

Discussion

Section V focused on practical design issues to consider when selecting patterns for structured illumination microscopy. Modeling the microscope OTF can be crucial to determining optimal pattern design. The model used here for example, has some known limitations as it is formed from an idealized symmetric PSF. A more complete model could take into account spherical aberrations of the objective lens and possibly effects of the specimen holder or covering member on the OTF. However, the results presented here demonstrate how to apply constraints from experimental parameters to design structured illumination patterns for microscopy tracking. Pattern selection can be further constrained by the requirement that the pattern be manufacturable with a high degree of precision to minimize tracking error. In one embodiment, larger scale patterns can be printed with photographic reduction to create the gradients required for sinusoidal patterns.

Photolithography can provide a method to produce smaller structures, usually characterized by sharp edges. Current research in microlithography extends the ability to print fine structures with gradients accurately with a reconfigurable, maskless approach [22].

Using a frequency selection method that maximizes the pattern SNR under experimental constraints, demonstrates lateral and axial tracking at or slightly above the noise floor. Different pattern designs and tracking strategies can work better under different scenarios. With for example, the simulated 4OX, 0.65NA objective lens, which has an Abbe resolution of 0.423 μm and depth of field of 2.34 μm. No specimen present: If a specimen is not present, pattern should contain for example, low frequencies. Lateral tracking with the MSC technique obtains an

RMS error of 1.1 nm (4.9 x 10^~3 pixel, 2.6 x 10^~3 Abbe units) when the pattern is in focus. Axial tracking can be accurate to within 0.50 μm (2.2 pixel, 1.2 Abbe units) when the pattern is less than 15 μm from the focal plane. Observing a specimen with sparse contrast: If observing a specimen having sparse contrast, pattern should contain for example, one low frequency to enable large displacements between frames. Several high frequencies that are above the dominant specimen frequencies should be present. If the pattern will remain less than 5 μm from the focal plane, lateral tracking with MSC obtains RMS errors better than 0.023 μm (0.10 pixel, 0.053 Abbe units). From 5 to 15 μm from the focal plane, WPC obtains RMS errors better than 0.046 μm (0.20 pixel, 0.11 Abbe units).

Observing a specimen with dense contrast: In this scenario, pattern should contain one low frequency for example, to enable large displacements between frames, and several high frequencies that are above the dominant specimen frequencies. Lateral tracking is accurate to within 0.11 μm (0.50 pixel, 0.27 Abbe units) if the pattern remains less than 5 μm from the focal plane. Above this, the SNR can drop too low to guarantee tracking accuracy. For axial tracking, the pattern should have a wide spread of frequencies, selected from local maxima in the SNR optimization. Tracking can be accurate to within 0.5 μm for axial distances of less than 15 μm, 6.4 times the depth of field.

Vl. Pattern removal In one embodiment, a goal of the system is to obtain high resolution microscopy mosaics. It can be crucial to remove the effect of pattern from the images. Otherwise, the structured illumination introduced by the pattern would interfere with the understanding of the structural details of the specimen. Pattern should be semitransparent so that the signal from the specimen behind the pattern can be recovered. A pattern that transmitted 50% of the light would result in a decrease of camera dynamic range by 1 bit. Because the pattern is a known transmission layer, its position and appearance can be detected by image analyzer, such as for example, software. Two approaches may be used.

Approach 1 : A Local Transmission Metric (LTM)

Consider the mean intensity image from all frames in a microscopy video. The mean intensity of a video with sufficient motion has low spatial variance, except in regions that contain stationary occlusions. The effect of motion on a mean intensity image is similar to the effect of a mean filter on a single image: the intensity from a single moving object is spread over many pixels in the mean image. Given sufficient uncorrelated motion, each point in the scene is imaged by all pixels within some neighborhood. The result is that the mean intensity image is approximately constant within local neighborhoods. For example, the predominant spatial variance in Figure 4.1e is due to stationary background objects. This observation leads to a method for estimating the stationary transmission map. Assume that for a sequence of N images the mean intensity,

is locally constant except where occluded by stationary background objects. If the occlusion regions are relatively small compared to the image plane, each occluded pixel will have within a small neighborhood a number of unoccluded pixels. The mean intensity in this neighborhood will be greater than the mean intensity at the occluded pixel. Comparing the mean intensity of an occluded pixel to that of its unoccluded neighbors determines the amplification necessary to bring the occluded pixel in line with its neighbors. In a sense, the occluded pixel appears as impulse noise, where the value that should have appeared at that pixel has been replaced by another value due to the occlusion. One common method used to detect impulse noise is to compare a pixel value to the median of neighboring pixel values [23]. For a neighborhood Q,

T_Ω(x) = Median[ϊ(Xi)] , {X; (X₁ GΩ} The local median of the equation above provides an estimate of how much light was transmitted to the pixels within a neighborhood. Knowing the expected size of occlusions in the image plane enables choosing a neighborhood size such that all neighborhoods have fewer than half of their pixels occluded. For this condition to be met for all neighborhoods, the occlusions present must be relatively small and sparsely distributed across the image plane. When such a neighborhood size can be chosen, however, an estimate of the stationary transmission map is:

This equation ensures that if a pixel has an intensity that is always lower relative to its neighbors, its intensity will be amplified. Similarly, if a pixel has an intensity that is always higher than its neighbors, its intensity will be diminished. Note that this latter case exposes a weakness in this approach applied to bright field microscopy: an occluding object would never amplify the intensity at a pixel, only diminish it. Approach 2: A Global Transmission Metric (GTM)

An alternative estimate of the stationary transmission map is using this method. For example, a stationary occlusion transmits a constant proportion of light at an image location. Foreground objects move over the occlusion throughout a video, so the intensity at that location changes over time. What remains constant, however, is the ratio of light transmitted at an occluded pixel compared to its neighbors. The transmission map estimate is constructed by forming a model of the stationary gradient of transmission at each pixel in the image plane. The transmission image formation: l(x; t) = L(X)T _c (X)T₅ (x; t), can contain two time-independent terms for example, L(x) being the spatially- varying illumination of the lamp and Tc(x) being the spatially-varying transmission from stationary components in the specimen. The nonuniform illumination can be calibrated and corrected using flat-fielding techniques. In the absence of at field calibration, however, these two terms are indistinguishable that is, both terms result in a stationary, nonuniform illumination at the image plane. The two terms, therefore, can be combined: l(x; t) = T_c (x)T_s (x; t).

As noted by Shizawa, the logarithm transforms this multiplicative imaging model into a linear imaging model [24]: log l(x; t) = log T_c (x) + log T_s (x; t). Taking the gradient of both sides of this equation yields:

V log l(x; t) = V log T_c (x) + V log T_s (x; t).

A benefit of transforming the imaging equation in this way is that moving components in the specimen affect the constant components only at the edges of the moving components. Said another way, if moving edges move sufficiently, there should be some estimator for which the contribution of the moving edges over all images is zero, for example:

E[V lOg T₈ (x; ti )] = 0, {i : i=1...N}. Applying such an estimator, E, to the equation above yields:

E[V Io₉ l(x; t; )] = E[V Io₉ T_c (x)], {i : i = 1...N}. The task remains to find a suitable estimator that can distinguish between elements from stationary occlusions and the moving specimen. Because the measurements of constant transmission at a single pixel can become corrupted with gross errors from moving transmission components, robust estimation suggests itself. In practice, the median can perform well in a number of cases, with the understanding that the median will provide a reliable gradient estimate at a pixel that is covered by moving foreground edges no more than half the time. Notice that the use of the gradient means that a large moving object may completely cover a region of the image plane in all frames, but as long as its edges are moving the constant transmission within this region can nevertheless be recovered. The logarithm of the constant transmission map is computed by integrating the gradient estimate over the image plane: logT_c(x) = JE[VlOgI(M_i)IdS ₁ {i : I = 1... N}, c where C denotes any path over the image plane from the origin to x. Converting back to intensity space can reveal the estimate of constant transmission,

T_c (S) = exp( JE[V log 1(Mi )]ds), {i : i = 1... N}.

C

The method of estimating the stationary transmission map from this equation is referred to as the gradient logarithm transmission (GLT) method and its variants due to estimator selection as GLT-mean and GLT-median.

Note that if the estimator, E, is linear, the estimation and gradient simplify the equation to:

T_c(x) - exp(E[logI(x,t_i)]), {i : i = 1...N]}

For the mean estimator, this equation reduces to a geometric mean of all frames. (As a result, the GLT-mean method does not actually involve a gradient computation in its implementation. No such simplification is possible for the median estimator.

Using the Approach- ImageTracker In one embodiment, ImageTracker is a software application written specifically for motion analysis in video microscopy. A chain of image filters can modify image sequence loaded into ImageTracker before further analysis. Filters include for example, intensity thresholding, Gaussian smoothing, gradient magnitude, and flat-fielding. ImageTracker can use the open-source Insight Toolkit (ITK) and Visualization Toolkit (VTK) to manage image loading, analysis, and display. This makes it easier to add more filters to future versions of ImageTracker.

As mentioned above, observations of in vitro specimens can for example, suffer from the specimen moving across the field of view, which makes analysis of motion within the specimen difficult. ImageTracker includes multi-resolution image registration that stabilizes image sequences with respect to this global motion. A user can select the coarsest and finest level of detail to use for the image registration, which enables targeting specific types of motion. For example, selecting the coarsest scale to smooth the majority of image features aligns large image structures. Selecting the finest scale to smooth small structures makes the alignment less sensitive to motion at that fine scale. After alignment, motion analysis of the fine structures can continue without concern for the global drift in the original video.

Long working distance lenses required for in vitro observation can result in large depths of field, which means specimen elements outside the target observation plane contribute to the video. When focusing through a non- moving layer it can be difficult to interpret the motion in a moving layer. ImageTracker can compute models of the non-moving layer by compiling statistics on the gradient of transmitted light over the entire image sequence and integrating the estimated non-moving transmission gradient. This model enables computing the image that would have been formed with only the moving layer present. This method can also correct for non-uniform illumination and debris on the image sensor.

Optical flow is a class of computer vision techniques for estimating the motion in an image sequence. Classic optical flow techniques can use either local image calculations or global energy. ImageTracker uses an implementation of an optical flow algorithm that combines local and global motion estimates. ImageTracker integrates frame-to-frame optical flow into net displacement fields. ImageTracker also visualizes the computed optical flow field along with the original image sequence. Two visualization techniques are provided.

Arrow glyphs can represent flow magnitude and direction at a sampling of image locations. Alternatively, a pseudo colored height map can represents flow magnitude over the entire image plane.

Discussion

This chapter presented two methods and software application for removing occlusions from microscopy videos that exhibit advantages over background subtraction. Each method and variation has strengths which can lead it to outperform the other in different scenarios. For example, in one embodiment the LTM method performs best on videos with small, bounded occlusions, but handles large occlusions and diffraction effects poorly. In one embodiment, the GLT methods can work equally well on any size occlusion and at least visually improves diffraction effects. In one embodiment, the GLT- median estimation can outperform the GLT-mean estimation when a moving scene object covers an image region for less than half of the video. As mentioned earlier, flat-fielding is a commonly-used technique in bright-field microscopy that counters the effect of non-uniformities in Kδhler illumination. Flat-fielding uses a specimen-free calibration image, or a set of images to create a map of how much light reaches each pixel in the image sensor. For example, a typical flat- field image is brightest at the center and falls off towards the edges. Flat-field correction can involve dividing specimen images by the flat-field calibration image. The transmission map computed in stationary occlusion removal is similar to the flat-field calibration image. In fact, because no distinction is made between light loss due to no uniform illumination and light loss due to fixed occlusions, stationary occlusion removal also performs flat-fielding. Given enough microscopy images in which no part of the specimen is stationary, stationary occlusion removal will recover the flat-field image for nonuniform illumination. Therefore, there can be two major advantages to stationary occlusion removal. First, stationary occlusion removal can be applied as a post-process in the absence of calibration data. Second, stationary occlusion removal recovers a specimen-specific calibration image that flat-fields non- uniformities from stationary parts of the specimen.

VII. Conclusion As mentioned earlier, one goal of methods, systems, and imaging media for microscopy tracking is to create a high-definition image mosaic of a large field of view, see Figure 38. In Figure 38, mosaic 300 can comprise a compilation of smaller images 302 which can for example, become pieced together using pattern data. The pattern data can be analyzed for lateral and axial positioning by using an image analyzer, such as computer imaging media using for example, frequency algorithms on the patterned specimen holder or covering member. The patterned specimen holder or covering member can allow for microscopy stage tracking, which can also eliminate the requirement of having costly motorized stages in microscopy. Once mosaic is formed, the pattern data can be removed.

The disclosure of each of the following references is hereby incorporated herein by reference in its entirety.

REFERENCES [1] P. Vuylsteke and A. Oosterlinck. Range image acquisition with a single binary-encoded light pattern. Pattern Analysis and Machine Intelligence, IEEE Transactions on, 12(2): 148-164, Feb 1990. 44, 70

[2] M. Noguchi and S. K. Nayar. Microscopic shape from focus using active illumination. In International Conference on Pattern Recognition, volume 1 , pages A: 147-152, October 1994. 44, 70

[3] Ramesh Raskar, Michael S. Brown, Ruigang Yang, Wei-Chao Chen, Greg Welch, Herman Towles, Brent Seales, and Henry Fuchs. Multi-projector displays using camera-based registration. In IEEE Visualization: Proceedings of the conference on Visualization '99, pages 161-168, Los Alamitos, CA, USA, 1999. IEEE Computer Society Press. 70 [4] S. S. Beauchemin and J. L. Barron. The computation of optical flow. ACM Computing Surveys, 27(3):433-466, 1995. 39, 40, 76

[5] J. P. Lewis. Fast normalized cross-correlation. Vision Interface, 1995. 77

[6] J. A. Nelder and R. Mead. A simplex method for function minimization. The Computer Journal, 7(4):308-313, 1965. 78

[7] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes in C: The Art of Scientific Computing. Cambridge University Press, 3rd edition, 2007. 55, 79, 93, 143, 147, 187

[8] H. Foroosh, J. B. Zerubia, and M. Berthod. Extension of phase correlation to subpixel registration. Image Processing, IEEE Transactions on, 11(3): 188-200, March 2002. 39, 80

[9] James W. Cooley and John W. Tukey. An algorithm for the machine calculation of complex Fourier series. Mathematics of Computation, 19(90):297-301 , 1965. 80

[10] Kenji Takita, Takafumi Aoki, Yoshifumi Sasaki, Tatsuo Higuchi, and Koji Kobayashi. High-accuracy sub pixel image registration based on phase-only correlation. IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, E86-A(8):1925-1934, 2003.40, 80, 131

[11 ] James D. Foley, Andries van Dam, Steven K. Feiner, and John F. Hughes. Computer Graphics: Principles and Practice. Addison-Wesley Publishing Company, Inc, 1997. 23, 26, 83, 130

[12] Frederic J. Harris. On the use of windows for harmonic analysis with the Discrete Fourier Transform. Proceedings of the IEEE, 66(1):51-83, 1978. 57, 84, 86, 95 [13] Emanuele Trucco and AIessandro Verri. Introductory Techniques for 3-D Computer Vision. Prentice Hail, 1st edition, 1998. 30, 38, 92, 99, 131

[14] Jim Pitman. Probability, chapter 3, pages 140-258. Springer, 1993. 99

[15] G. Healey and R. Kondepudy. Radiometric ccd camera calibration and noise estimation. IEEE PAMI, 16(3):267-276, 1994. 26, 27, 30, 49, 58, 99, 113

[16] Ted PeIIa. Gilder TEM grids, commercial internet web site, available at http://www.tedpella.com/grids html/gilder. htm., March 2009. 111

[17] J. A. Hartigan and M. A. Wong. Algorithm as 136: A k-means clustering algorithm. Applied Statistics, 28(1): 100-108, 1979. 113

[18] W. Y. Hsu, W. F. Paul Poon, and Y. N. Sun. Automatic seamless mosaicing of microscope images: enhancing appearance with color degradation compensation and wavelet-based encoding. Journal of Microscopy, 231 (3):408- 418, 2008. 124

[19] ED Salmon, SL Shaw, JC Waters, CM Waterman-Storer, PS Maddox, E Yeh, and K. Bloom. A high-resolution multimode digital microscope system. Methods in Cell Biology, 56:185-215, 2003. 6, 125

[20] LA Mangiamele and SS Burmeister. Acoustically evoked immediate early gene expression in the pallium of female tungara frogs. Brain, Behavior and Evolution, 72:239-250, 2008. 151

[21] Theodore George Rochow and Paul Arthur Tucker. Introduction to Microscopy by Means of Light, Electrons, X Rays, or Acoustics. Plenum Press, 2nd edition, 1994. 155 [22] Melanie V. Kessels, Marwa El Bouz, Robin Pagan, and Kevin Heggarty. Versatile stepper based maskless microlithography using a liquid crystal display for direct write of binary and multilevel microstructures. Journal of Micro/Nanolithography, MEMS and MOEMS, 6(3):033002, 2007. 158

[23] R. Garnett, T. Huegerich, C. Chui, and Wenjie He. A universal noise removal algorithm with an impulse detector. Image Processing, IEEE Transactions on, 14(11):1747-1754, Nov. 2005. 53

[24] Masahiko Shizawa and Kenji Mase. Simultaneous multiple optical estimation. ICPR, 1 :274-278, 1990. 42, 50, 54

Claims

CLAIMS What is claimed is:

1. A system for microscopy tracking, the system comprising: a specimen holder or covering member for holding or covering a specimen being imaged using a microscope; and a pattern formed in or on the specimen holder or covering member, the pattern being structured for tracking, through at least one of: frequency, phase, and defocus analysis, microscope stage position for images of the specimen obtained at different stage positions of the microscope.

2. The system of claim 1 , wherein the specimen holder or covering member comprises a specimen slide.

3. The system of claim 1 , wherein the specimen holder or covering member comprises a specimen slide cover.

4. The system of claim 1 comprising an image analyzer for automatically tracking the stage positions using the pattern.

5. The system of claim 4, wherein the image analyzer is configured to track the stage positions in lateral directions using spatial frequencies in the pattern.

6. The system of claim 4, wherein the image analyzer is configured to track the stage positions in an axial direction using defocus analysis.

7. The system of claim 4, wherein the image analyzer is configured to use the pattern to form a mosaic image from the images.

8. The system of claim 4, wherein the image analyzer is configured to remove the pattern from the images.

9. The system of claim 1 , wherein the pattern is structured to encode absolute lateral stage position via frequency and phase of the pattern structure.

10. A method for microscopy tracking, the method comprising: obtaining images of a specimen obtained at different stage positions of a microscope, where the images include images of a pattern located in or on a specimen holder or covering member, where the pattern is structured to allow tracking of the different stage positions through at least one of frequency, phase, and defocus analysis; and tracking the stage position using the at least one of frequency, phase, and defocus analysis for each of the images.

11. The method of claim 10, wherein the specimen holder or covering member comprises a specimen slide.

12. The method of claim 10, wherein the specimen holder or covering member comprises a specimen slide cover.

13. The method of claim 10, wherein tracking the stage position includes tracking the stage position in lateral directions.

14. The method of claim 13 wherein tracking the stage position in lateral directions includes tracking the stage position using spatial frequencies in the pattern.

15. The method of claim 10, wherein tracking the stage position includes tracking the stage position in an axial direction.

16. The method of claim 15, wherein tracking the stage position in the axial direction includes tracking the stage position using defocus analysis.

17. The method of claim 10, comprising using the pattern to form a mosaic image from the images.

18. The method of claim 10, wherein the pattern is structured to encode absolute stage position through structures of different frequencies in the pattern and wherein determining the stage position includes determining an absolute stage position by using phase analysis of the encoded structure.

19. The method of claim 10, comprising removing the pattern from the images using an image analyzer.

20. A computer readable medium having stored thereon executable instructions that when executed by the processor of a computer control the computer to perform steps comprising: obtaining images of a specimen obtained at different stage positions of the microscope, where the images include a pattern located in or on a specimen holder or covering member, where the pattern is structured to allow tracking of microscope stage position through at least one of frequency, phase, and defocus analysis; and analyzing the images using the at least one of frequency, phase, and defocus analysis to determine stage position for each of the images.