968

J. Opt. Soc. Am. A / Vol. 31, No. 5 / May 2014

Trattner et al.

Image formation of thick three-dimensional objects in differential-interference-contrast microscopy Sigal Trattner,1,5 Eugene Kashdan,2,4,* Micha Feigin,3,4 and Nir Sochen4 1

Department of Medicine, Division of Cardiology, Columbia University Medical Center, New York 10032, USA 2 School of Mathematical Sciences, University College Dublin, Belfield, Dublin 4, Ireland 3 MIT Media Lab, Massachusetts Institute of Technology, 75 Amherst St., Cambridge, Massachusetts 02139, USA 4 Department of Applied Mathematics, Tel Aviv University, Ramat Aviv, Tel Aviv 69978, Israel 5 Department of BioMedical Engineering, Tel Aviv University, Ramat Aviv, Tel Aviv 69978, Israel *Corresponding author: [email protected] Received November 4, 2013; revised January 27, 2014; accepted February 26, 2014; posted March 19, 2014 (Doc. ID 200319); published April 9, 2014 The differential-interference-contrast (DIC) microscope is of widespread use in life sciences as it enables noninvasive visualization of transparent objects. The goal of this work is to model the image formation process of thick three-dimensional objects in DIC microscopy. The model is based on the principles of electromagnetic wave propagation and scattering. It simulates light propagation through the components of the DIC microscope to the image plane using a combined geometrical and physical optics approach and replicates the DIC image of the illuminated object. The model is evaluated by comparing simulated images of three-dimensional spherical objects with the recorded images of polystyrene microspheres. Our computer simulations confirm that the model captures the major DIC image characteristics of the simulated object, and it is sensitive to the defocusing effects. © 2014 Optical Society of America OCIS codes: (110.2990) Image formation theory; (180.0180) Microscopy; (290.0290) Scattering. http://dx.doi.org/10.1364/JOSAA.31.000968

1. INTRODUCTION The differential-interference-contrast (DIC) microscope is a powerful noninvasive interferometer for visualizing biological specimen. It is a light microscope supplemented with the DIC module, which translates phase variations into intensity images. The DIC microscope enables visualization of transparent objects that otherwise cannot be seen in the standard light microscope, unless stained. A general goal of microscopy is to extract quantitative information from the recorded images, specifically the shape and the measures of the visualized object. The shape reconstruction of the object is obtained by analyzing sets of two-dimensional images made using different focus settings (so-called “optical sectioned images”). Usually the reconstruction is preceded by the removal of image degrading effects introduced by the optical system. In computational optical sectioning microscopy (COSM), the shape reconstruction is usually done by the deconvolution methods. The idea behind these methods is in modeling image intensity as a convolution of the object’s intensity transmittance with a computed point-spread function (PSF). (PSF is a response of an imaging modality to a point source or an object point.) A linear function of image intensity enables us to obtain a high-quality restored image in modalities such as fluorescence and bright-field microscopy. However, in DIC microscopy, the task of shape reconstruction is challenging due to nonlinear connection between physical properties of object and image intensity as well as due to directional sensitivity yielded by the optical system. The problem becomes even more difficult due to cumulative build-up of scattering mass if the observed specimen are thick 1084-7529/14/050968-13$15.00/0

(more than 20 μm in diameter) [1–4]. Thus, developing a methodology for computer simulation of the image formation process in DIC microscopy, specific to thick objects, has become a problem worth special attention. Solution of this problem and its subsequent experimental verification are the goals of this manuscript. A structure of DIC microscope is shown in Fig. 1. A DIC microscope is a visible light microscope combined with a polarizer–analyzer pair and two Wollaston prisms. These additional components add interferometric functionality to the microscope, which in turn enables the visualization of transparent (phase) objects [4–6]. A number of DIC image formation models are proposed and discussed in the literature [7–11]. The majority of these models are limited to thin objects. In most cases, the first-order Born approximation [12] is used for modeling of light scattering. However, as has been shown in [13–15], the first-order (linear) Born approximation and similarly defined Rytov approximation are inapplicable to the image formation modeling of thick biological specimen. For instance, they cannot be used in the study of human embryo cells, which could reach 100 μm in diameter. The high-order Born and Rytov approximations that can be used to model more complicated scattering patterns experience severe convergence issues. In particular, detailed analysis provided in [17] shows that the Born series diverges when the refractive index of the object is different by more than 20% from the environment. The high-order Rytov approximation experiences very similar issues. In addition, the computer approximation of the Rytov series is made especially difficult as the scattered potential does not have a finite support as it does in the Born series [16]. © 2014 Optical Society of America

Trattner et al.

Vol. 31, No. 5 / May 2014 / J. Opt. Soc. Am. A

Detector Analyzer 2nd Wollaston Prism Objective Object

Condenser 1st Wollaston Prism Polarizer

Light Source Fig. 1. Diagram of principal components of DIC microscope.

The product of convolutions method is often used for image formation modeling in phase contrast microscopy [18]. In [19], the authors applied it to DIC microscopy. In this method, a three-dimensional object is modeled as a stack of parallel planes transverse to the optical axis. This approach excludes contributions from backscattering and relies on the discretization along the optical axis, which limits its applicability to simulate the scattering exhibited by thick objects. The geometrical optics approach to modeling DIC image formation can be attractive due to its simplicity and low computational complexity [9]. However, the axial focus distortion—a major effect in the image formation process of thick objects—cannot be accurately modeled by geometric ray-tracing methods [20]. In general, it can be stated that geometrical optics can sometimes provide a good approximation to wave optics problems if the light is incoherent in both time and space [21]. A different approach to numerical modeling of DIC image formation uses vectorial Maxwell’s equations to compute the scattering of light [25]. However, this approach was computationally tested on small (less than 1 μm in diameter) particles. As a result, the applicability of this approach to visualization of thick objects remains unproven. The finite difference time domain (FDTD) algorithms became the most common way of solving Maxwell’s equations to compute light scattering by general geometries [28]. Yet, when the simulated objects are thick, using FDTD becomes cumbersome due to the size of the problem, despite existing commercial and in-house software written for this purpose. The implementation of these methods requires resolving issues that are beyond the scope of our paper. In particular, FDTD is based on rectangular grid and requires a special approach to deal effectively with curved interfaces. The complexity increases due to the need to sample the entire computational domain rather than only the object. Finally,

969

the absorbing boundary conditions should be formulated and implemented to avoid nonphysical wave reflections [22,23,25]. Image formation analysis of thick objects entails a number of computational challenges, which have to be answered in order to obtain meaningful results from the experimentalist point of view. Among the major problems encountered are sampling (limited by the minimal resolution requirements from above and exponentially raising computational cost from below), approximation of highly oscillatory integrals, round-off error accumulation due to the processing of large volumes of data, and convergence issues. The image formation model validity for thick objects remains questionable until neither of these problems is met and addressed properly. To the best of our knowledge, no approach to DIC image formation was validated on thick objects (with diameter d ≫ 10 μm), representing a very important sector of DIC microscopy applications. In this manuscript, we present an image formation model that simulates electromagnetic wave scattering and propagation through the DIC microscope components, and validate it against actual images of thick transparent objects. Our approach to DIC image formation modeling is designed for objects of general shape and has no size limitations. However, in this work, we refer to a homogeneous (uniform refractive index) sphere as our object of interest in order to use the closed form of solution to Maxwell’s equations given by the Lorentz–Mie theory and run proper model validation through comparison with the experiments. The feasibility study of the model is a necessary requirement to its extension to “real-life” problems, and it prerequisites application of other computational techniques such as the FDTD to simulation of light scattering by general three-dimensional shapes. Tshe ultimate goal of our DIC image formation model will be its application to formulation and solution of the inverse problem dedicated to extraction of volumetric information from a series of DIC images. The manuscript is organized as follows: in Section 2 we describe our DIC image formation model corresponding to the different stages of light propagation through the microscope components. In Section 3, we explain the experimental part of our work—recorded measurements of the polystyrene microspheres used for model verification. Section 4 includes a brief overview of the computational approach we use for simulation of the DIC image formation and our considerations behind the choice of the numerical algorithms employed for this purpose (the details of the computational algorithm and its implementation are put into the corresponding appendices). Results of numerical simulations compared both qualitatively and quantitatively with the recorded DIC images are displayed and analyzed in Section 5. The manuscript is concluded by the discussion in Section 6.

2. PROCESS OF DIC IMAGE FORMATION A. Physical Background The light route in the DIC microscope begins with a single light source, usually from the filament of a tungsten–halogen lamp. The light is polarized by the polarizer and then splits into two light waves by the first Wollaston prism. These waves are mutually coherent, orthogonally polarized, laterally displaced by an extremely small distance (called the shear), and phase-shifted relative to each other (a shift called the bias

970

J. Opt. Soc. Am. A / Vol. 31, No. 5 / May 2014

Trattner et al.

module, which consists of the polarizer and the first Wollaston prism and emanates as two light waves that enter the condenser (Fig. 3). We consider an incident light wave from a common source as a monochromatic, coherent, electromagnetic field, Ei
r (the superscript i stands for incident). Its amplitude is defined as E0 , and it could be expressed as

retardation). After passing through the condenser, the light waves travel through the closely adjacent regions of the simulated object. Each wave is phase-shifted depending on its optical path through the object. The light waves are assembled and focused by the objective lens before entering the second Wollaston prism to be reunited into a single wave. The next placed analyzer cancels the orthogonality between the components of the united waves and enables their interference. The amplitude, and thus the intensity of the resulting wave, is a function of the phase difference (gradient) between two waves. Therefore, the difference in the optical path length of the neighboring points in the object produces a variance in the intensity of the image, namely, the contrast. It enables the visualization of outlines in the transparent object by a CCD camera or by the human eye [2]. The intensity variation image reflects phase differences and makes the DIC microscope suitable for visualization of transparent objects, in particular, living specimen and cells. The phase gradients yield a shadowing artifact (the relief-like structure that DIC images are known for). The physical formation of the DIC image is broken into a sequence of elements. Following the physical path of the light waves through the DIC microscope components, we refer to the five key transitions. For each transition we identify the major physical process, model it, and compute the corresponding electric fields over the relevant plane, perpendicular to the optical z axis. The model consists of the following stages corresponding to the DIC microscope elements:

Ei
r  E0 exp
ikˆez · r − iωt;

where r 
x; y; z are the Cartesian spatial coordinates with the z axis lying along the light propagation, eˆ z . Other parameters include the angular wave number of the incident field k and the fixed angular frequency ω. Time variation of the electric field is omitted further in the text, yet assumed as exp
−iωt (the field is considered time-harmonic). In actual physical reality, perfect coherence does not exist, and real sources of light are never monochromatic. The choice of monochromatic light is motivated by analytical purposes and mathematical convenience. And yet, it can be justified. Subject to the restriction that an arbitrary electromagnetic wave can be Fourier analyzed into a superposition of plane monochromatic waves, a solution to the interaction of light and object can be obtained by the superposition of fundamental solutions. This linearity is attributed to Maxwell’s equations and the boundary conditions [12,24]. It is also claimed that the basic monochromatic theory is adequate enough to describe the action of interference modalities [12]. Preza et al. [10] extended the scalar DIC model to partially coherent illumination and showed agreement of experimental results with the full coherence assumption. The incident light is unpolarized when entering the polarizer. It becomes polarized according to the direction of the polarizer transmission axis, φ, and the light wave is depicted as

1. Light transition through the first module of the DIC — the polarizer and the first Wollaston prism — and the condenser; separation of the single light wave into two separate waves with different phases (Section 2.B). 2. Scattering of two light waves by the simulated object (Section 2.C). 3. Assembling and focusing of light by the objective lens (Section 2.D). 4. Propagation of light between the objective and the image plane (Section 2.E). 5. Transition of light waves through the second module of the DIC — the analyzer and the second Wollaston prism; intensity calculation of the DIC image (also included in Section 2.E).

Ei
r  E 0 exp
ikˆez · rˆeφ ;

1 E
x e ex ; z · rˆ 1
r  E 0 exp
ikˆ 2 1 e−z · rˆey : E
y 1
r  E 0 exp
ikˆ 2

B. Transition of Light through the First DIC Module The route of a light wave in the DIC microscope originates at the light source. The light wave passes through the first DIC 1st Wollaston Prism

Polarizer

E (i) eˆx

Condenser

Objective E

E1(x) (y) 1

E

(3)

Note that the waves are angularly split by a small angle, the shear angle. After passing through the condenser the waves

Object

E (i)

(2)

where eˆ φ is the unit vector in the φ direction. The polarized light then strikes the first Wollaston prism. This birefringent device splits the light wave to receive two mutually coherent, orthogonal, plane-polarized light waves. Without loss of generality, we take x, y to be the polarization directions. The polarization planes of two waves make an angle of 45° with the incident wave polarization plane, eˆ φ , as illustrated in Fig. 4, and can be expressed as

A diagram of the suggested model is charted in Fig. 2. The Cartesian coordinate system is localized such that the z axis coincides with the optical axis of the microscope. The axis origin is set at the center of the simulated-visualized object.

Light Source

(1)

(x) 2

E2(y)

y

(x) 3

E

E3(y)

2nd Wollaston Analyzer Prism E4(x)

E

E5 Optical Axis

(y) 4

x z

Detector

Image plane

Fig. 2. Diagram of the DIC image formation process with the field notations at key planes.

Trattner et al. Light Source

E (i)

Vol. 31, No. 5 / May 2014 / J. Opt. Soc. Am. A 1st Wollaston Prism

Polarizer (i)

E eˆx

B⃗  μH⃗

Condenser (x) E1

Object δr

(y) E1

Optical Axis

Fig. 3. Transition of light through the first DIC module.

emanate as plane waves [26] and travel in parallel paths. The angular split is transformed into extremely small spatial distance, the shear distance δr. The light wave fields can be now written as 1 E
x ez · rˆex ; 1
r  E 0 exp
ikˆ 2 1 ez ·
r  δrˆey : E
y 1
r  E 0 exp
i
kˆ 2

(4)

A bias retardation phase shift, ϕbias , can be added to one of the light waves by adjusting the Wollaston prism. The two light waves incident upon the object can then be generally summarized as 1 E
x ez · rˆex ; 1
r  E 0 exp
ikˆ 2 1 ez ·
r  δr  ϕbias ˆey : E
y 1
r  E 0 exp
i
kˆ 2

(5)

We refer to E
x 1
r as representing the parallel polarization and to E
y
r as the perpendicular polarization in regard to 1 the x–z plane [27]. From this point and on we continue and depict the model referring to a single wave. We take E
x 1
r as the parallel-polarized field for describing the propagation of light in the system. The same rules will be valid for E
y 1
r, with the adequate adaptations for spatial and phase differences. These rules will be described later on, when both waves will be considered again at the interference stage (Section 2.E). C. Scattering of Light by the Object Following light interaction with the object we obtain the scattered electromagnetic field, which is described by Maxwell’s equations for electric E and magnetic H fields in the isotropic medium: ∂B⃗  ∇ × E⃗  0; ∂t

∂D⃗ ⃗ − ∇ × H⃗  −J; ∂t

(6)

where J is an electric current density. Electric and magnetic fields are connected with the corresponding fluxes (D and B) through the constitutive relations:

Fig. 4. Polarization planes of the split waves with respect to the incident field direction.

⃗ D⃗  εE;

971

(7)

where ε is the dielectric permittivity and μ is the magnetic permeability. In our model we simulate pure phase objects that are visualized; therefore the light absorption is negligible and the dielectric coefficient is considered a function of space only ε  ε
⃗r . As stated in Section 1, our main interest is to simulate all stages of the DIC image formation process for thick objects and validate our approach experimentally. Therefore, we consider spherical objects as an initial example for our purpose due to the closed formulation of the time-harmonic electromagnetic field scattered by the spherical objects given by the Lorenz–Mie theory (Mie series solution). The solution of time-harmonic Maxwell’s equations by Mie series assumes that the object is composed of a homogeneous, isotropic, and optically linear material irradiated by an infinitely extending plane wave. Our implementation of Mie series solution is based on the algorithms discussed in [27,29]. The details of the mathematical description of the Lorenz–Mie theory can be found in the literature, in particular, in [12,27,30,31]. We outline the Mie series approach using the notation of Bohren and Huffman [30] in Appendix A. Since an arbitrarily polarized plane wave can be expressed as a superposition of two orthogonal polarization states, the scattering problem should be solved twice (for a given direction of propagation) [30]. In the DIC microscopy, the wave separation is convenient as the incident light wave naturally splits into two orthogonal polarized waves, passing through the first Wollaston prism. After passing through the object, the light field can be described using two expressions for the incident light and the scattered light, denoted Ei2 , Es2 , respectively. The expressions for both fields are derived in Appendix A using Lorenz–Mie theory. These fields are received upon the lens, and their superposition forms the total field E2 —the result of the propagation of the parallel-polarized field E
x 1
r through the object. D. Field Transformation through the Objective A use of infinity-corrected lenses in light microscopy has grown rapidly in recent years due to multiple industrial and research applications [32]. The infinity-corrected objective is composed of an objective chosen by the microscope user and a relay lens, called the tube lens, which is fixed. The relay lens is located within the body tube, between the objective and the intermediate image (which we refer to as the image plane). A thin lens approximation provides a good basis for modeling relay lens transformation. The thin lens approximation idealizes the process of light passage and is commonly used for microscope modeling [9,33]. A more rigorous model for infinity-corrected objectives was suggested by Meinhart and Wereley [34]. And yet, they have shown that under moderate magnification, the diffraction-limited spot size of this model differs very slightly from the diffraction-limited spot size of the thin lens model, implying that thin lens modeling can be exploited for most applications. We implement the thin lens model of the aberration-free ideal lens as it is presented by Goodman [24]. The thin lens model is based on the paraxial approximation. Small angles are assumed; namely, light from a point (xo , yo ) on the object is close to the optical axis in the sense that xo ∕zo and yo ∕zo are small enough.

972

J. Opt. Soc. Am. A / Vol. 31, No. 5 / May 2014

Trattner et al.

The infinity-corrected lens modeled as a thin lens is shown in Fig. 5. The field, E2 , received upon the lens as discussed in Section 2.C, is passed through the lens by the thin lens transformation function T l [24] E3
u; v  T l E2
u; v;

(8)

where u and v represent the coordinates of the thin lens model aperture and E3 stands for the vector fields received on the exit from the lens. The lens center is located at the point where the optical axis crosses the lens. A phase transformation introduced by the lens depends on the variable phase delay, which is a function of the distance from the lens center. The incident and the scattered fields on the exit from the lens are denoted Ei3 and Es3 , respectively. The transformation function T l of the lens is expressed as
T l
u; v  exp

 −ik
u2  v2  P
u; v; 2f

1 1 1   ; f zo zi

(10)

where zo represents the distance between the object and the lens and zi is the distance between the lens and the image plane. P
u; v is the pupil function of the circular lens, representing the limits of the finite extent of the lens aperture: P
u; v 

1; u2  v2 ≤ R2 ; 0; otherwise;

Objective

E 2(x)

E4
ξ; η 

E3(x)

eikzi iλzi

  ik
u − ξ2 
v − η2  dudv; E3
u; v exp 2zi −∞ (12)

ZZ

∞

where ξ, η are the polar coordinates of the detector. The finite limits of the aperture are already incorporated in E3 . We use the vector Kirchhoff approximation [35]; i.e., each Cartesian component of the field is individually subjected to a scalar diffraction calculation. Thus, for each Cartesian component of the field E4 , the diffraction integral is solved separately. They solutions are then recombined at the end, to form the corresponding vector field. The propagated diffracted incident and scattered fields yield the total field for each of the polarized components:

(11)

where R is the lens aperture radius, which we define as the tube lens radius. In order to describe correctly the compound objective, we choose the thin lens parameters in accordance with the focal lengths of the infinity-corrected lens components. The distance from the object to the lens, zo , is defined as the focal length of the objective, f lens . The distance from the lens to the image plane, zi , is equal to the focal length of the tube lens, f tube . In general, light waves passing through the lens are changing both their phase and polarization. However, modern DIC specialized lenses are designed to preserve the polarization by removing stress and strain from the glasses of the numerous lens elements [32]. Therefore we assume that the polarization

Object

E. Propagation of Light from the Lens to the Image Plane Following the exit from the microscope objective lens, the light waves fall on the image plane. We ascribe the image plane of the DIC microscope to the second Wollaston prism, the analyzer, and the detector, as shown in Fig. 6. The light waves pass these components sequentially, in that order. The propagation of the total field E3 from the lens to the image plane is computed using the two-dimensional Fresnel diffraction integral [24] over the lens aperture:

(9)

where f is the focal length of the thin lens calculated from the thin lens equation:



directions of two wave fields passing the lens are unaffected by the transmission through the lens and only the phase is changed.


xi
xs E
x 4
ξ; η  E 4
ξ; η  E 4
ξ; η;

(13)


yi
ys E
y 4
ξ; η  E 4
ξ; η  E 4
ξ; η;

(14)

where E
x 4
ξ; η is the parallel-polarized diffracted field and E
y 4
ξ; η is the perpendicular-polarized diffracted field. These wave fields then enter the second Wollaston prism. The second prism is similar to the first prism but it is oriented in a manner that enables combining two waves to a single one. The shear distance between the two waves that was previously introduced by the first prism (Section 2.B) is canceled by the second prism. The recombined wave still has two perpendicular polarization directions. The wave then passes through the analyzer (a second polarizer) that projects the wave on a single polarization plane. The interference is revealed in this plane. Optimally, the transmission plane of the analyzer should be normal to the polarizer transmission plane for pure phase objects. As shown in Fig. 7, the analyzer transmission plane makes angles of 45° and 135° with the polarization planes 2nd Wollaston Analyzer Prism

Objective

Optical Axis z0

y x

zi

E 2(y)

(y) E3

Image Plane

z Fig. 5. Diagram of the thin lens approximation: f is the focal length, zo is the object-to-lens distance, and zi is the lens-to-image plane distance.

Object

(x) 3

E

E4(x)

E3(y)

E4(y)

Detector

E5 Optical Axis

zi

Image Plane

Fig. 6. Light diffraction and transition through the second DIC module.

Trattner et al.

Vol. 31, No. 5 / May 2014 / J. Opt. Soc. Am. A Polarizer Transmission plane

E4(x)

By computing field intensity we complete the formulation of the DIC image formation algorithm.

3. EXPERIMENTS

E4(y) 135° 45°

Analyzer Transmission plane

Fig. 7. Configuration of analyzer with respect to vibration planes of the incoming waves.

of the parallel-polarized field E
x and the perpendicular4 polarized field E
y 4 , respectively. According to Fig. 7, the field interference is equal to the sum of the projected fields on the analyzer plane for each Cartesian component:
y E5
ξ; η 
E
x eφ⊥ ; 4 cos
π∕4  E 4 cos
3π∕4ˆ

(15)

where φ⊥  φ  π∕2, φ is the transmission plane of the polarizer (see Fig. 4) and φ⊥ is the transmission plane of the analyzer. Note that the frequency dependence was deferred until this stage, since we have dealt with the monochromatic light. Yet, since the intensity is the squared amplitude of the field, it is necessary to incorporate the contribution of all frequencies, before its calculation. On this stage we restore the notation of Section 2.B and add the dependency on the frequency ω to the field E5 . At this point the contribution of each frequency yields a fundamental solution, E5
ξ; η; ω. Using the superposition principle and according to the light source profile p
ω the resulting field can be expressed as Z E5
ξ; η 

∞ 0

p
ωE5
ξ; η; ωdω:

(16)

In computations, the integral is approximated by the finite sum over the contributing frequencies of the visible spectrum shown in Fig. 8. The field intensity I can then be derived as 2 2 2
y
z I
ξ; η  jEj2  jE
x 5
ξ; ηj  jE 5
ξ; ηj  jE 5
ξ; ηj : (17)

100 90 80

% Transmission

973

70 60 50 40 30 20 10 0 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700

Wavelength [nm]

Fig. 8. Transmission curve of the band-pass filter.

A. Microscope and CCD Setup The microscope used in our experiments is a transmitted light Zeiss AxioImager M1 microscope with the DIC module. The objective used for acquiring the images is a multi-immersion Plan Neufloar 25X∕0.8 DIC with a pupil diameter of 10 mm. The condenser used is an Achromatic-aplanatic universal condenser 0.9H D of Zeiss. The light source in the microscope is HAL 100 halogen illuminator, 12V. A band-pass filter was localized in front of the light source to limit the spectral bandwidth for practical modeling. This filter has a peak at 0.549 μm and a full width at half-maximum (FWHM) equal to approximately 78 nm. The catalog number of the filter has not been provided by the manufacturer; therefore the transmission curve of the filter was simulated in a spectrophotometer (Cary 5000) and is displayed in Fig. 8. In this microscope, the analyzer and polarizer are fixed to a right, a setting that is optimal for imaging phase objects as it simplifies the interpretation of the DIC images [6] (see also Section 2.E). Bias retardation phase could be introduced by sliding the Wollaston prism. In the microscope used for our experiments, the bias retardation phase had no calibrated setting. Therefore, the images of darkest background were considered to have no bias retardation phase added, i.e., ϕbias  0 radians. The brightest background images were considered to have a maximal added bias retardation phase equal to ϕbias  π radians [6,10]. Bias retardation phase between
−π; 0 radians could be gained by sliding the Wollastion prism (using a nonscaled screw). The model was examined on thick transparent polystyrene microspheres of diameter 82 μm, with refractive index of n  1.59–1.60 (at wavelength λ  0.589 μm). The microspheres were immersed in deionized water [indexed 7732-18-5] with refractive index n0  1.33 (Polybead microspheres of Polysciences, Inc., catalog No. 07315). The images were recorded by a camera AxioCam MRm from Zeiss with 12 bit depth and a well size of 6.45 μm × 6.45 μm. Table 1 summarizes the major settings of the DIC image acquisition experiment. The microspheres were put on a standard microscope slide of approximately 1 mm thickness and covered with a cover glass of (0.17  0.005 μm (Superior Marienfeld cover glass No. 1.5H). A small chamber was sealed to avoid the microsphere squeezing. High vacuum grease was used to make the chamber’s walls. Since the microspheres are embedded in the deionized water, water was also used as an immersion medium for the lens to minimize spherical aberrations. The optical components (light source, filter, objective, and condenser) were aligned to provide coherent Köhler illumination. (Köhler illumination is a method for providing object illumination that is uniformly bright and free of glare; for more information, see [12,38].) Coherent illumination is achieved when the condenser is fully closed so its numerical aperture (NA) equals zero [39]. In practice, the aperture of the condenser could not be fully closed, but was brought to a minimum. A closed condenser aperture enhances the image contrast but reduces the NA of the condenser and thus the resolution. The microscopic resolution (distance between two

974

J. Opt. Soc. Am. A / Vol. 31, No. 5 / May 2014

Table 1. Summary of the Experimental Setup of the DIC Microscope and the Microspheres Viewed Microscope Components Microscope Zeiss Axioimager M1 Condenser lens 0.8-NA Band-pass filter Peak at 0.549 μm, FWHM 78 nm DIC Objective Multi-Immersion Type Plan Neufluar Imm Korr DIC Numerical aperture 0.8-NA Magnification 25X Pupil diameter 10 mm Immersion used Water Camera Parameters CCD camera Zeiss Axiocam 12 bit 12 V CCD well 6.45 μm × 6.45 μm Microsphere Parameters Microsphere diameter 82 μm Material Polysterene n  ∼1.59–1.60 Immersion material Deionized water n0  ∼1.33

neighboring points), in our case, is roughly equal to the single wavelength (less than 1 μm). Such resolution is sufficient for our purposes due to the large size of the simulated object [40]. B. Image Acquisition and Correction The microscope and camera were controlled by the AxioVision microscope software package of Zeiss. Automatically regulated focus enabled the acquisition of optical sectioned images in different depths along exact intervals of the z axis. The spacing interval between the optical sections (the focal planes along the z axis) was set to 0.5 μm. Illumination controlling software enabled us to use the maximal dynamic range of the camera without saturation. In order to avoid saturation of the camera CCD wells, the illumination exposure time was regulated per image set. Various “parasitic” effects might occur during the image acquisition process, and digital image restoration should be performed to eliminate the corresponding artifacts [39]. Removal of the systematic noise and the artifact reduction from the DIC images is addressed in Section 4.B. To simplify the comparison between the modeled and the recorded DIC images and to reduce the nonuniformity of the illumination, the experimental images were subtracted from the background. An “empty” image without the microspheres served as a background on each level of focus. The background image had similar illumination conditions and system parameters to the corresponding image, which included the microspheres. Thus, any potential nonuniformities of light, varying sensitivity of the detector, or artifacts could be compensated for [39].

Trattner et al.

Table 2. Summary of Image Formation Model Parameters, Corresponding to the Experimental Setup Parameter Spherical object diameter Refractive index of spherical object Refractive index of medium Wavelength range Lens–detector distance Focal length of thin lens Thin lens aperture diameter

Value 82 μm n  1.60 n0  1.33 0.501–0.580 μm 164.5 mm 6.6 mm 3 mm

medium surrounding them (water). The modeling of the objective compound lens according to the thin lens equation (10) requires two parameters, zo and zi (Fig. 5). The distance between the thin lens and the image plane (zi ) was set as 164.5 mm. The objective focal lens (f obj ), reported by the manufacturer as 6.6 mm, was used as the distance between the object and the lens (zo ). The thin lens aperture diameter parameter was set as 3 mm; the pupil diameter was reported by Zeiss as 10 mm. To reduce computational complexity we used a smaller aperture diameter than reported by the manufacturer. The simulated images are computed as the weighted superposition of the fundamental solutions for different incident wavelengths, according to Eq. (16). The wavelength interval was determined by the band-pass filter and the weight attributed to each wavelength. Due to practical considerations only the wavelengths with significant effect on the image were taken into account. Computationally this means that 80 wavelengths in the range of 0.501–0.580 μm are sufficient to form the image. The parameter values used for computation of the model are summarized in Table 2. B. Data Adaptation In order to compare recorded images to the model-based simulations, we incorporate several physical effects into the simulated images that were not taken into account in the “ideal” mathematical model derived in Section 2. In particular, we account for the microscope and the CCD noise inherent in the

4. NUMERICAL SIMULATION OF DIC IMAGES A. Model Parameters The image formation was simulated using the model described in Section 2 on the example of a dielectric microsphere with parameters that corresponded to the experimental setup. The light scattering by the sphere was modeled using the Mie series solution as it is discussed in Appendix A. The simulated microsphere diameter was set as 82 μm, according to the measurements on the recorded images. The refractive indices were set as n  1.60 for the spheres and n0  1.33 for the

Fig. 9. Effect of adding blur and noise to the simulated image: diagonal profile of the modeled “ideal” image (dashed blue line), and diagonal profile of the modeled image after convolution with low-pass filter with added Gaussian white noise, N (0, 0.0003) (solid red line).

Trattner et al.

Vol. 31, No. 5 / May 2014 / J. Opt. Soc. Am. A

system. We made two model adjustments that refer to the CCD readout noise and the diffraction of the compound objective. Zero-mean Gaussian noise can be seen as fluctuating small oscillations riding on the signal envelope. The white Gaussian noise is associated with the CCD readout noise that originates from the readout amplifier of the CCD, and that is considered as a significant noise component [39]. The quality of DIC images is also degraded by the blurring effects due to the diffraction of the compound lens (occurring in fact in any microscopical modality [10]). Since the objective lens of the microscope is modeled using the thin lens model, the diffraction of the system is missing. We convolved a low-pass kernel with the simulated image to compensate for this effect. The low-pass kernel is smoothing the profile

975

by reducing the sharpness of the image. One could replace this approach by truncating the field in the pupil planes of the condenser and objective; however, the result is very similar to that of convolution with a low-pass filter, and we used it for convenience. The convolution is simple, fast, and straightforward from the numerical point of view. Figure 9 shows the effect of both the blur and the white noise addition on the diagonal profile of the simulated image.

5. RESULTS A. Images with no Bias Retardation Phase Recorded optical sectioned DIC images of the polystyrene microsphere were compared to the modeled images. As

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Fig. 10. Recorded (left) and simulated (center) DIC images of 82 μm sphere with no bias: (a)–(c) in focus, (d)–(f) 20 μm above focus, (g)–(i) 40 μm above focus, and (j)–(l) 20 μm below focus. The condenser aperture is closed to maximal position; the images are processed for the system noise; the right column shows the comparison of the diagonal profile of the recorded (dashed blue line) and simulated (solid red line) images.

976

J. Opt. Soc. Am. A / Vol. 31, No. 5 / May 2014

Trattner et al.

less steep along with the distance, as can be observed in the summarizing graph of recorded diagonal profiles in Fig. 11(a). The colors red, black, and blue represent 10, 20, and 40 μm object distances above the focus level, respectively. On the other hand, the corresponding profiles of the simulated images do not show widening, and changing slopes are not evident at first sight due to the ringing effect [Fig. 11(b)]. Yet, the trend observed in Fig. 11(a) could be clearly seen in the simulated images after disregarding the ringing effect as shown for the 20 μm above focus image [illustrated by the black dashed line in Fig. 11(b)]. B. Images with Bias Retardation Phase A different setting of the microscope takes into consideration an added bias retardation phase. The bias retardation phase was set to π∕4 radians for modeled image calculations and subsequent comparison to the recorded images. A series of recorded DIC images of the microsphere, with a bias addition, is shown in Fig. 12 for different focus levels. The recorded images are shown with the corresponding simulated images with a bias addition. For each pair of images a diagonal profile across the shear direction is presented. The blue line represents the profile of the recorded image, and the red line represents the profile of the simulated image. Recorded DIC images are characterized by one main lobe with high intensity and a dark, low-intensity, lobe-like part across the shear direction. The addition of bias retardation phase emphasizes the optical bias-relief structure for which the DIC images are known. The intensity of the brighter lobe decreases towards the center of the microsphere. The image is asymmetric with respect to the shear direction, and a strip with brighter intensity is located at the dark lobe exterior part. As the image is recorded away from the in-focus location, the edges become more blurry and expanded. The background of the simulated image is brighter than the background of the corresponding recorded image with no bias added. Comparison between simulated and recorded DIC images shows that the model captures major DIC image characteristics. Both the bright lobe and the dark lobe are captured, as well as the bright strip. The contrast and the intensity ratio of the simulated images are in good qualitative agreement with the recorded images. The defocusing effect of the external part of two lobes is also observed in the simulated images. However, the ring-like effect is present in the simulated images.

1

1

0.9

0.9

0.8

0.8

Normalized intensity

Normalized intensity

described in Section 4.B, both recorded and simulated images were processed to reduce the system noise. The details on model implementation and on numerical algorithms used for computations are described in Appendix B. We call the image in focus if the equator of the microsphere is located at the lens focal plane. In the case of the recorded images, the in-focus image is selected visually and it serves as the anchor reference for all other optical sectioned images (located from both sides of the focal plane). Figure 10 displays results for different focus levels with no bias retardation. Each figure shows the normalized intensity of a simulated image versus the corresponding recorded image. Diagonal profiles are drawn for each couple of images from both the recorded and the simulated images. The diagonal profiles are taken along the shear direction of the DIC images (45° to the image horizontal axis). The blue lines represent profiles of the recorded image, and the red lines represent profiles of the simulated image. Recorded images of the DIC microsphere are characterized by two main lobes across the shear direction with high intensity at the edges and decreasing intensity towards the center of the circle-like image of the microsphere. The image is highly symmetric with respect to the shear direction. When the image is recorded away from the in-focus location, the edges become more blurry and expanded at their external part. The comparison between simulated and recorded DIC images shows that the model captures the main DIC characteristics of two sheared lobes as well as the gradual decrease towards the center. The normalized intensity of the simulated images shows an agreement in contrast and intensity ratio with the recorded images. The defocusing effect of the twolobe external part is also observed in the simulated images. Yet, some of the simulated images suffer from the ring-like effect. It is seen both in the inner part of the microsphere as well as at the edges. These could be explained by the discrete nature of the model, system aberrations that are not included in the model, and other effects related to the “ideal”nature of the model, which are discussed in detail in Section 6. Simulated images of the microspheres placed above infocus level have very similar size compared to the recorded microsphere. On the other hand, the simulated images of the microspheres located below the focal plane look bigger in size than the corresponding experimental measurements. Both distinct sides of the recorded microsphere image located above the focal plane become wider, and the slope becomes

0.7 0.6 0.5 0.4 0.3

0.7 0.6 0.5 0.4 0.3

0.2

0.2

0.1

0.1

0

10 20 30 40 50 60 70 80 90 100 110 120 130 140

0

10 20 30 40 50 60 70 80 90 100 110 120 130 140

Length [micrometers]

Length [micrometers]

(a)

(b)

Fig. 11. Diagonal profiles of (a) recorded DIC images and (b) simulated DIC images at different distances above focus; no bias retardation. Solid red, dotted black, and angled dashed blue lines represent 10, 20, and 40 μm above focus level, respectively. The dashed black line represents the profile of the 20 μm above focus simulated image following removal of the ringing effect.

Trattner et al.

Vol. 31, No. 5 / May 2014 / J. Opt. Soc. Am. A

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

977

Fig. 12. Recorded (left) and simulated (center) DIC images of 82 μm sphere with bias retardation: (a)–(c) in focus, (d)–(f) 20 μm above focus, (g)–(i) 40 μm above focus, and (j)–(i) 20 μm below focus. The condenser aperture is closed to maximal position; the images are processed for the system noise. The right column shows a comparison of the diagonal profile of the recorded (dashed blue line) and simulated (solid red line) images.

6. SUMMARY AND DISCUSSION In this paper we have presented an image formation model for transparent thick spherical objects observed in a DIC microscope. The light propagation through the microscope components and its scattering by the object were modeled using a vectorial approach. The scattering pattern was computed using the Mie series solution of time-harmonic Maxwell’s equations. The model was evaluated by a comparison of simulated images with recorded DIC images of 80 μm diameter polystyrene microspheres. We have shown that the suggested model captures the major characteristics of a DIC image of three-dimensional thick objects. In particular, the shear effect of the DIC microscope is visible on the simulated images. Our DIC image formation

model is the first vectorial model that was successfully validated against experimental images of three-dimensional objects with diameter much bigger than a wavelength. Simulating image formation of thick objects was computationally challenging. Numerical implementation of the mathematical model entailed high oscillations of the wave function due to the ratio between the wavelength and the object size. Computing Mie series coefficients required large numbers of terms due to slow convergence. To speed up the computations we used a graphic processing unit (GPU), courtesy of NVIDIA. In this stage of the research, we omitted quantification of some of the simulation results and concentrated on a rather qualitative analysis. In particular, we did not quantify the

978

J. Opt. Soc. Am. A / Vol. 31, No. 5 / May 2014

Trattner et al.

changes in intensity, which could be normalized without compromising the quality of the model. Mathematical modeling of complex physical processes requires a number of assumptions and simplifications and involves computational shortcuts. Therefore the full “recreation” of the DIC image on the computer screen is quite difficult. For instance, one of the issues that requires further investigation is related to missing gray background on the modeled bias retardation phase images. Still, we have succeeded to capture the major features of the DIC image of thick transparent objects. We should notice that the ringing effect (so-called “Gibbs ringing”) that raises due to the truncation of Mie series and affects the quality of the simulated image should be reduced with minimal compromise of the image resolution. Two strategies we consider worth future investigation for this purpose are (i) post-processing based on morphological image analysis and subsequent removal of periodic (linear in the polar system of coordinates) features associated with the rings [41] and (ii) incorporation of the windowing function (i.e., the Hann filter) in the Fresnel integral (12) for ring removal [35]. Our testing of the second strategy yielded visibly reduced blurring effects, and an increase in the image contrast, but negligible effects on the rings. On the other hand, the resulting image was very sensitive to small shifts in the weight parameter of the filter, and further investigation is required [42]. Other features that would be addressed in the next version of the model will include aberrations by the objectives and polarization changes that occur as the light passes through the object and the microscope components. In real life, specimen are not totally transparent and possess some extent of light absorption that will be also taken into account. This work has proved the feasibility of the vectorial approach to DIC image formation modeling, in particular, for analysis of thick transparent objects that is required in multiple practical applications. It is a first major step towards a more general model that would represent the morphology of objects of higher complexity, i.e., biological cells. Our method sets a ground for future work on the corresponding inverse problem of extracting volumetric information from the DIC images, for instance, by use of a combination of the image formation model with the nonlinear optimization methods and subsequent error minimization.

APPENDIX A: MIE SERIES SOLUTION We outline the Mie series approach for parallel polarization. We review the expressions suitable for a sphere within a nonabsorbing medium only. Parallel (or x) polarization refers to a polarization parallel to the x–z plane. We consider the x-polarized plane wave incident on a homogeneous, isotropic sphere. Pursuing commonplace notations in the literature and the notation of Eq. (5), it is denoted i E
x ez · rˆex ; 1
r  E
kr  E 0 expikˆ

(A1)

where E 0 replaces
1∕2E 0 with respect to the notation in Eq. (5) for future convenience. We convert to spherical polar coordinates r, θ, and ϕ. The expansions of the incident electric field, the internal field, and the scattered field into the spherical harmonics are derived as

Ei
kr 

∞ X


1 E n M
1 o1n − iN e1n ;

n1

Eint
mkr 

∞ X


1 En cn M
1 o1n − idn N e1n ;

n1

Es
kr 

∞ X


3 E n ian N
3 e1n − bn M o1n ;


A2

n1

where an , bn , cn , and dn are unknown scattered and internal field expansion coefficients, respectively (for parallel polarization). The ratio between the refractive indices of the sphere, n, and the embedding medium, n0 , is denoted as m, and En is given as E n  in


2n  1 E : n
n  1 0

(A3)


1;3 The expressions for the spherical harmonics M
1;3 o1n and N e1n are
1;3 M
1;3 eθ o1n
r; θ; φ  cos φπ n
cos θzn
ρˆ

− sin φτn
cos θz
1;3 eφ ; n
ρˆ N
1;3 e1n
r; θ; φ  cos φn
n  1 sin θπ n
cos θ  cos φτn
cos θ − sin φπ n
cos θ

zn
1;3
ρ eˆ r ρ

0 ρz
1;3 n
ρ eˆ θ ρ

0 ρz
1;3 n
ρ eˆ θ ; ρ

(A4)

where πn 

P 1n
cos θ ; sin θ

τn 

dP 1n
cos θ ; dθ

(A5)

and P 1n
cos θ are the associated Legendre functions of the first kind of degree n and order 1; eˆ r , eˆ θ , and eˆ ϕ are the unit vectors in spherical coordinates; ρ  kr for both the incident field and the scattered field; and ρ  mkr for the internal field. Superscripts appended to M and N denote the spherical Bessel function zn . zn
1 denotes j n
ρ, which is defined as  0.5 π j n
ρ  J n1∕2
ρ; (A6) 2ρ where J n1∕2 is the Bessel function of the first kind. zn
3 denotes h
1 n (called the spherical Hankel function), where h
1 n
ρ  j n
ρ  iyn
ρ;

yn
ρ 

 0.5 π Y n1∕2
ρ; 2ρ

(A7)

(A8)

and Y n1∕2 is the Bessel function of the second kind. At the sphere–medium interface, a boundary condition is imposed due to the discontinuity at the boundary points, as a result forcing continuity of the tangential components of the electromagnetic field. In order to satisfy the boundary conditions, the electric fields must fulfill

Trattner et al.

Vol. 31, No. 5 / May 2014 / J. Opt. Soc. Am. A


Ei  Es − E int  × eˆ r  0;

(A9)

with the radial unit vector eˆ r . In physical terms, this boundary condition is a sufficient condition for energy conservation across that boundary [30]. The coefficients an , bn , cn , and dn are retrieved from the boundary conditions and computed using the Riccati–Bessel functions, ψ n
ρ  ρj n
ρ;

ξn
ρ  ρh
1 n
ρ;

(A10)

by substituting r  a, where a is the sphere radius, in Eq. (A2) and in the following expressions: mψ n
mkrψ 0n
kr − ψ n
krψ 0n
mkr ; mψ n
mkrξ0n
kr − ξn
krψ 0n
mkr ψ
mkrψ 0n
kr − mψ n
krψ 0n
mkr bn  n ; ψ n
mkrξ0n
kr − mξn
krψ 0n
mkr mξn
mkrψ 0n
mkr − mξ0n
mkrψ n
mkr cn  ; ψ n
mkrξ0n
kr − mξn
krψ 0n
mkr mξ0n
mkrψ n
mkr − mξn
mkrψ 0n
mkr dn  : mψ n
mkrξ0n
kr − ξn
krψ 0n
mkr

an 


A11

After passing through the object, the light field can be expressed in two expressions for the incident light and the scattered light, denoted Ei2 , Es2 , respectively. The fields are received upon the lens, and their superposition makes the total field E2 , which is the result of the propagation of E
x 1
r, the parallel-polarized field, through the object.

APPENDIX B: MODEL IMPLEMENTATION Implementation of the DIC image formation model has required us to deal with various computational challenges. The large size of the studied object with respect to the wavelength made the implemented numerical algorithms both complex and time consuming. The convergence properties of the algorithm and minimization of round-off errors arising during the computations required special attention. Computational time had to be manageable in order to introduce the model into practical use. A substantial amount of computation was required for the modeling of two major optical processes: the Mie series solution and the Fresnel diffraction integral. The infinite Mie series, in Eq. (A2), is truncated to the sum over a finite number of nc terms. However, the number of terms required for convergence can still be quite large, and round-off error accumulation can become a problem [30]. Error accumulation might yield incorrect results, in particular, distortion of the final image. The subtle trade-off between convergence of the algorithm and the round-off error was constantly examined. An analysis of the convergence behavior of the Mie series led to the following empirical rule for estimation of the number of terms [27,29]: nc  x  4.05x1∕3  2;

(B1)

where x  kr is the size parameter (k is the wave number and r is the radius of the sphere). Yet, this rule has some exceptions, e.g., the case of the spherical resonance mentioned by Barber and Hill [27].

979

Computing Mie series coefficients involves spherical Bessel functions and their derivatives. Due to the large size of the object and consequent potential problems of memory overflow, instability, and inaccuracy, these functions were approximated using recursion. The corresponding discussions and algorithms, including the approach we have implemented, can be found in the literature [27,29,30]. The Fresnel integral in expression (12) is highly oscillatory. To avoid loss of accuracy due to dense equidistant sampling of the aperture, we applied a multivariate Gaussian quadrature technique to the numerical integration [36,37] instead of the commonly used Newton–Cotes formulae. The first version of the program was written in MATLAB and was run on the Intel Core 2 Quad Q9550 CPU at 2.83GHz. The simulation time of a single wavelength solution exceeded 100 h. GPUs lead the way with the current trend of highperformance computing hardware. Current GPUs are massively parallel single instruction multiple data (SIMD) machines. A parallel version of the code that employed CUDA-based multithreading was executed on the NVIDIA Tesla S1070 GPU server. The running time for a single wave solution of the optimized for GPU code was approximately 80 min, yielding an image of 491 × 491 pixels [43]. Mie calculation is a pointwise operation and thus transfers well to the GPU. The algorithm has two parts: a preliminary calculation of the coefficients for the Hankel functions (Bessel functions of the first and second kind) and the associated Legendre functions, and the main part of the Mie computation. For the first part, we need to deal with proper memory management and mixed precision computations to reduce memory overhead. For the second part, we need to deal with numerical stability issues in order to allow for single precision computations enforced by GPU. Tables 3 and 4 illustrate the CPU and wall time required for calculations of the Mie series solution and Fresnel integral [44]. To maintain speed, while keeping accuracy in check, mixed accuracy methods are often required. Even when computations using floating point numbers seem to be problematic, enough accuracy can be achieved with care. There are two optimizations that helped us achieve enough stability for using floating point numbers: (i) avoiding adding a large number to a small number and (ii) calculating modulo the period to avoid computing trigonometric functions of large argument, which yields loss of significant digits.

Table 3. Run Times for Mie Series Computations, in hh:mm:ss

CPU time Wall time

Intel Core 2 Quad Q9550

Tesla S1070 GPU

5:59:12 5:53:46

0:03:05 0:02:40

Table 4. Run Times for Fresnel Integral Computations, in hh:mm:ss

CPU time Wall time

Intel Core 2 Quad Q9550

Tesla S1070 GPU

166:42:15 63:10:39

1:22:10 1:21:34

980

J. Opt. Soc. Am. A / Vol. 31, No. 5 / May 2014

Trattner et al.

ACKNOWLEDGMENTS The authors are grateful to Dr. Ben Ovryn for the fruitful discussion about the image formation modeling in light microscopy and to the anonymous reviewers for their suggestions, which helped us to improve the manuscript significantly.

20. 21. 22.

REFERENCES 1. 2. 3. 4. 5. 6. 7. 8.

9. 10.

11.

12. 13. 14.

15.

16. 17. 18. 19.

D. Agard and J. Sedat, “Three-dimensional architecture of a polytene nucleus,” Nature 302, 676–681 (1983). S. Bradbury and P. Evennett, Contrast Techniques in Light Microscopy. Microscopy Handbooks 34 (Bios Scientific, 1996). J. Padawer, “The Nomarski interference-contrast microscope. An experimental basis for the image interpretation,” J. R. Microsc. Soc. 88, 305–349 (1967). M. Pluta, Advanced Light Microscopy, Vol. 2 (Elsevier Science, 1988). W. Lang, “Nomarski differential interference contrast microscopy. I. Fundamentals and experimental designs,” Zeiss Information 70, 114–120 (1968). W. Lang, “Nomarski differential interference contrast microscopy. II. Formation of the interference image,” Zeiss Information 71, 12–16 (1969). M. Arnison, K. Larkin, C. Sheppard, N. Smith, and C. Cogswell, “Linear phase imaging using differential interference contrast microscopy,” J. Microsc. 214, 7–12 (2004). N. Axelrod, A. Radko, A. Lewis, and N. Ben-Yosef, “Topographic profiling and refractive-index analysis by use of differential interference contrast with bright-field intensity and atomic force imaging,” Appl. Opt. 43, 2272–2284 (2004). F. Kagalwala, F. Lanni, and T. Kanade, “Computational model of DIC microscopy: from observations to measurements,” Technical report CMU-R1 TR (Carnegie Mellon University, 2000). C. Preza, D. Snyder, and J. Conchello, “Theoretical development and experimental evaluation of imaging models for differentialinterference-contrast microscopy,” J. Opt. Soc. Am. A 16, 2185–2199 (1999). E. Van-Munster, L. Van-Vliet, and J. Aten, “Reconstruction of optical pathlength distributions from images obtained by a wide-field differential interference contrast microscope,” J. Microsc. 188, 149–157 (1997). M. Born and E. Wolf, Principles of Optics, 4th ed. (Cambridge University, 1999). S. Trattner, M. Feigin, H. Greenspan, and N. Sochen, “Validity criterion for the Born approximation convergence in microscopy imaging,” J. Opt. Soc. Am. A 26, 1147–1156 (2009). S. Trattner, M. Feigin, H. Greenspan, and N. Sochen, “Can Born approximate the unborn? A new validity criterion for the Born approximation in microscopic imaging,” in Mathematical Methods in Biomedical Image Analysis (MMBIA) Workshop, in conjunction with ICCV’07, Rio de Janeiro, Brazil (2007). S. Trattner, E. Kashdan, M. Feigin, M. Greenspan, C.-F. Westin, and N. Sochen, “DIC microscopic imaging of living cell and error analysis of Born approximation,” in Proceedings of 3rd Workshop on Microscopic Image Analysis with Applications in Biology, in conjunction with MICCAI’08 (2008). A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (IEEE, 1988). M. Slaney, “Imaging with diffraction tomography,” Ph.D. thesis (Purdue University, 1985). H. Sierra, C. A. DiMarzio, and D. H. Brooks, “3D effects in DIC images of extended objects,” Proc. SPIE 7184, 71840D (2009). H. Sierra, C. A. DiMarzio, and D. H. Brooks, “Modeling phase microscopy of transparent three-dimensional objects: a

23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42.

43.

44.

product-of-convolutions approach,” J. Opt. Soc. Am. A 26, 1268–1276 (2009). J. J. Stamnes, Waves in Focal Regions (Adam Hilger, 1986). P. Török, S. J. Hewlett, and P. Varga, “The role of specimeninduced spherical aberration in confocal microscopy,” J. Microsc. 188, 158–172 (1997). E. Kashdan and E. Turkel, “A high order accurate method for the frequency domain Maxwell’s equations across interfaces,” J. Sci. Comput. 27, 75–95 (2006). E. Kashdan and E. Turkel, “High order accurate modelling of electromagnetic wave propagation across media: grid conforming bodies,” J. Comput. Phys. 218, 816–835 (2006). J. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996). R. Munro and P. Török, “Vectorial, high numerical aperture study of Nomarski’s differential interference contrast microscope,” Opt. Express 13, 6833–6847 (2005). M. Shribak and S. Inoué, “Orientation-independent differential interference contrast microscopy,” Appl. Opt. 45, 460–469 (2006). P. Barber and S. Hill, Light Scattering by Particles: Computational Methods (World Scientific, 1990). A. Taflove and C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, 2005). W. Wiscombe, “Improved Mie scattering algorithms,” Appl. Opt. 19, 1505–1509 (1980). C. Bohren and D. Huffman, Absorption and Scattering of Light by Small Particles (Wiley, 1983). M. Mishchenko, L. Travis, and A. Lacis, Scattering, Absorption and Emission of Light by Small Particles (NASA Goddard Institute for Space Studies, 2006). MicroscopyU, http://www.microscopyu.com/. S. Schaub, D. Alexander, and J. Barton, “Theoretical model of the laser imaging of small aerosols: applications to aerosol sizing,” Appl. Opt. 30, 4777–4784 (1991). C. D. Meinhart and S. T. Wereley, “The theory of diffractionlimited resolution in microparticle image velocimetry,” Meas. Sci. Technol. 14, 1047–1053 (2003). B. Ovryn and S. Izen, “Imaging of transparent spheres through a planar interface using a high-numerical-aperture optical microscope,” J. Opt. Soc. Am. A 17, 1202–1213 (2000). P. David and P. Rabinowitz, Methods of Numerical Integration (Academic, 1975). W. Gautschi, Orthogonal Polynomials: Computation and Approximation (Oxford University, 2004). O. M. Primer, http://micro.magnet.fsu.edu/primer/. S. Inoué and K. Spring, Video Microscopy: The Fundamentals, 2nd ed. (Plenum, 1997). H. H. Hopkins and P. M. Barham, “The influence of the condenser on microscopic resolution,” Proc. Phys. Soc. B 63, 737–744 (1950). J. Sijbers and A. Postnov, “Reduction of ring artifacts in high resolution micro-CT reconstructions,” Phys. Med. Biol. 49, N247–N253 (2004). S. Trattner, E. Kashdan, H. Greenspan, and N. Sochen, “Human embryo under the DIC microscope—vectorial approach to the electromagnetic scattering simulation,” in Proceedings of 8th International Conference on Spectral and High-Order Accurate Methods (ICOSAHOM), Trondheim, Norway (2009). S. Trattner, M. Feigin, E. Kashdan, and N. Sochen, “GPU accelerated electromagnetic scattering and diffraction in 3D microscopic image formation,” in Proceedings of the 3rd Workshop on GPUs for Computer Vision, Barcelona, Spain (2011). M. Feigin, “Computational methods in image analysis,” Ph.D. thesis (Tel Aviv University, 2012).