Multiple angle digital holography for the shape measurement of the unpainted tympanic membrane.

The shape of the tympanic membrane (TM) plays an important role in sound transmission through the ear for hearing. Previously we developed a high-speed holographic system employing a tunable wavelength laser for rapid TM shape measurement. However, the tunable laser illumination was not sufficient to measure the shape of the unpainted TM due to the semi-transparency of the TM and short exposure time of the camera. This paper presents a new multiple angle illumination technique that allows us to use a higher power single wavelength laser to perform shape measurements on the unpainted TM. Accuracy of the new method is demonstrated by a measure of a step gauge provided by the National Institute of Standards and Technology. We successfully applied the new shape measurement method on a fresh postmortem human TM without any paint.


Introduction
Contactless and whole-field measurement of the surface shape with high resolution is required in various fields of science and technology. Digital holographic techniques have a three-dimensional imaging capability and, thus, becomes an interesting alternative for shape measurements of diffusively reflecting surfaces. Surface shape can be computed from the difference of the reconstructed phase recorded before and after the change of either incident angle of illumination wave [1][2][3][4][5][6][7][8] or wavelength of the used light [9][10][11][12][13][14][15][16][17][18]. A specific example is the measurement of the tympanic membrane (TM) shape for hearing research. The TM couples sound energy from the environment to the inner ear [19,20]. The transmission of energy through the middle ear is affected by the TM's complex geometry (shape and thickness), spatially varied mechanical properties and microstructure [21][22][23][24][25][26][27][28][29][30][31]. The shape information together with sound induced displacements of the TM are essential to understanding mechanics of the TM. Moreover, the shape of the TM also has a diagnostic value for middle ear pathologies, e.g. bulging of the TM is a sign of middle ear effusion. Measurements of the TM shape have been carried out utilizing various techniques, including Moiré interferometry [32], optical coherence tomography [33], ultrasound [34], digital image processing [35] and digital holography (DH) [8]. Our group has developed a high-speed multiple wavelength DH employing a tunable laser to measure the TM shape in cadaveric human ears [36][37][38][39]. However, the tunable laser illumination employed in our system was not sufficient to measure the shape of the unpainted TM due to the semi-transparency of the TM and short exposure time of the camera. Painting the TM is not practical for in vivo shape measurement. To overcome this limitation, an alternative approach is to use a higher power single wavelength laser to resolve the shape with the two-angle DH method [8]. However, the conical-shaped nonhomogeneous and semi-transparent properties of the unpainted TM can result in strong optical phase field decorrelation between two even slightly changed illumination angles, which limits accuracy and robustness of the two-angle method. In this paper, we introduce a new multiple angle illumination technique that allows us to use a single wavelength (532nm) higher power laser to perform robust and accurate shape measurements on the unpainted TM surface. Accuracy of the new method is estimated by a measure of a step gauge provided by National Institute of Standards and Technology. We successfully applied the new shape measurement method on a fresh postmortem human TM without any paint. Moreover, the developed method is not limited to TM shape measurement. A similar method has been used to measure shape and hence acoustic properties of a curved glass shell with attached piezoelectric macro fiber composite (MFC) actuators [40], which can be further adapted to describe the precise placement of elements within optical systems.

Multiple angle digital holography
Digital holography (DH) is a technique that records and reconstructs an optical wave field. The optical wave is described by its amplitude and phase (polarization state is omitted). However, only the amplitude information can be retrieved by optical detectors as the measured intensity I. Due to the lack of phase-sensitive detectors, it is necessary to encode the phase information into the measurable intensity data. Such an encoding is based on superposition of an object wave carrying information about the measured quantity with a reference wave. The interference pattern H is called a digital hologram and can be described by the interference formula [41] where A and B stand for the additive and multiplicative terms, respectively. The interference phase P represents phase differences between the object wave and the reference wave. The hologram H is stored in a computer memory as an array of numbers. In digital holography, a digital hologram H represents the wavefield in the hologram (detector) plane. This wavefield can be numerically propagated using diffractive optics theory [42] resulting in a complex field U that represents the wavefield at the image plane, where the reconstructed image is focused. Phase P and intensity I information can be retrieved as [41] P = tan −1 [Im{U}/Re{U}] and I = |U| 2 In case of recording the wavefield at the image plane by digital holography, a lens focusing the object directly on a digital camera is employed. Therefore, the hologram and the image planes coincide and no wavefield propagation is needed. This approach is interchangeable with speckle interferometry. Then the complex field U is obtained by e.g. a simple four step phase shifting technique [42], where four π/2 phase shifted digital holograms H 1 , H 2 , H 3 , H 4 are captured. The complex field U can be computed using the formula where i = √ −1. The phase shifting technique allows the in-line interferometer configuration and thus the use of the full pixel count in forming the holographic image.
Let us consider the geometry of the sample in Fig. 1. The optical path difference between the reference beam path SR 1 R 2 D and the object beam path SOD is related to the interference phase The interference phase includes the shape information represented as variation in the z coordinate. Samples to be measured are diffusively reflecting surfaces and hence digital holograms H 1 , H 2 , H 3 , H 4 as well as the interference phase P exhibit speckle noise. The speckle noise appears due to random interference of the scattered wave and its appearance is almost independent of the object characteristics, but strongly depends on the optical properties of the viewing system. Therefore, shape z cannot be retrieved from a single hologram. Small variations in illumination angle dθ, however, result in a phase difference: That difference is used for measurement of the surface shape. Speckle patterns of the two phase fields (before and after the angle change) are supposed to be correlated and suppressed as they are subtracted from each other. Due to the harmonic nature of light, the phase difference is wrapped within the range of 2π radians. The measured phase difference dP includes a linear term that carries no information about the sample shape. On the other hand, the shape z (after the removal of the linear term dP TLT ) can be retrieved as: Assuming a sensitivity factor dk = (2π/λ)(sin θdθ) and substituting it into Eq. (7), the phase difference simplifies to: From Eq. (8) one can see that larger variations in dθ of the object illuminating beam increase the sensitivity factor dk, which brings better measurement accuracy. On the other hand, with higher sensitivity, surface contours having distance ∆z = λ/(sin θdθ) create denser fringe patterns that can be difficult to spatially unwrap. Moreover, the speckle correlation is lower and phase difference dP can be noisy (especially for non-homogeneous and semi-transparent surfaces like TM). Gradients of the measured surface also limit the maximal applicable dθ and sensitivity factor dk, respectively. Theoretically, a maximal phase jump between any two adjacent pixels of dP must be smaller than half of the contour distance, i.e. λ/(2 sin θdθ), in order to avoid an ambiguity. Practically, the slope should be even smaller due to speckle noise. To gain maximal accuracy while avoiding the ambiguity or slope issues, a multiple angle technique is proposed. As the illumination angle θ sweeps, a set of n=1,2,. . . N interference phases P n (θ) is computed. The sequence P n (θ) is processed in a differential way (using two consecutive interference phases) so that: where dθ m =θ n+1 -θ n and m is an integer m=1,2,. . . N-1. Considering only z-coordinates (after removal of the term dP TLT ), Eq. (11) becomes: depth z_max, i.e.: ∆z > z max . Even if the phase differences d P m must be spatially unwrapped, the unwrapping is much more robust and unlikely to fail due to the very low phase range of d P m (only a few fringes are expected in the phase difference maps). The next step proceeds with a cumulative sum of the phase difference sequence d P m : resulting in a stack of phase difference maps. The unambiguous phase difference sequence dP m has the same meaning as dP in Eq. (9), however, the phase range is no longer wrapped into the range of 2π radians. The object shape z can be computed as a slope of the unambiguous phase difference with respect to the vector of the sensitivity factor, see Eq. (9):

Experimental arrangement
The experimental setup is illustrated in Fig. 2. A laser beam was generated by Nd:YAG laser (LAS) of wavelength λ = 532.8 nm and maximal output power of 50 mW. The laser beam was split using a beam splitter (BS) into reference (RB) and object (OB) beams. After reflection at mirror M4, the object beam propagates through a microscope objective (MO1) converting the beam into a spherical wave. The spherical wave reflects from mirror M5 placed on a rotational stage (RS) and illuminates a sample. The sample is located L 2 =445 mm from the vertex point (VP) of the mirror M5 and the illumination angle is θ=32°. These values were determined using yardstick measurement. Scattered light from the object is collected by the camera lens (CL) and imaged on a high-speed camera (CAM) sensor. The high-speed camera-Photron Fastcam SA-Z with a pixel pitch of 20 µm operates at 67.2 kHz (512×512 pixels). The employment of a high-speed camera minimizes low-frequency natural noises in a living sample, such as motions due to respiration, heartbeat, muscle tremor, etc. and thus makes quantitative measurements on live subjects possible. The reference beam is phase shifted by movements of mirror M3 that is attached to a piezoelectric transducer (PZ) and diverged using another microscope objective (MO2). High speed cameras have usually lower resolution. To keep the full spatial resolution, the in-line interferometer configuration together with the phase shifting technique is applied. For each angle, a sequence of at least four phase shifted holograms must be captured. Due to the use of the high speed camera, such a sequence is captured in less than 1 ms. The interference pattern (digital hologram) of both superposed beams (RB + OB) recombined by wedged beam splitter (WBS) is captured by the CAM. During the measurement, the rotation stage (RS) with mirror M5 is continuously rotating from its initial position (mirror position MP1) β=0°to its final position (MP2) at approximately β≈0.8°. The rotation takes about 1s. At several time instants (further denoted as mirror positions) during the RS movement, a sequence of phase shifted holograms is captured and processed using Eqs. (3) and (2) in order to compute phase difference maps. The sequence of the phase difference maps is processed using Eqs. (11)- (14) to retrieve the surface shape. If the angle of rotation was bigger, higher sensitivity of the method would be achieved, but at the same time, the light beam displaces out of the field of view. This could be improved by modification of the experimental arrangement. Assuming maximal angle span allowed by hardware, the mirror is tilted at an angle β≈0.8°and the central ray of the beam is reflected at angle α=2β≈1.6°. It is important to note that due to the intended application for in-vivo measurement inside the narrow ear canal (see Fig. 13), the illumination beam is a spherical diverging wave; not a plane wave. As the mirror rotates, virtually, the source point (S1) of the diverging beam shifts (to S2), see dashed area in Fig. 2. Assuming the distance from the vertex point (VP) to the source point (S1): L 1 =|VP S1|=22 mm and the angle α≈1.6°, the virtual shift is a=0.6 mm. Multiple angle contouring method presented in this paper measures optical phase difference between the two optical wavefields before (position MP1) and after (position MP2) the source point shift. Unlike contouring using plane waves, the interference fringes are curved creating the arcs of hyperbolae. However; we can assume that the source points are localized symmetrically on the axis of dθ and that the source points lie far from the sample i.e. L 1 +L 2 a. In this case it simply follows that the fringes are rectilinear and vertical as in the case of two plane waves contouring at angle: The angle span of the mirror rotation β≈0.8°, as well as other parameters of L 1 and L 2 leading to computation of the illumination beam angle change dθ, is estimated by trigonometric computation from the yardstick measurement. However; due to the small rotation angle, such measurements have large uncertainty. Moreover, for the presented temporal data processing, values of the angle steps within the full angle span must be known. Therefore, accurate values of rotation are determined using standard flat surface sample as described in the next section.

Flat surface sample calibration
Before TM shape measurement, the system was calibrated by measuring a flat sample fabricated from black anodized aluminum. Flatness of the surface within field of view of the camera is better than 5 microns and thus does not significantly contribute to the uncertainty budget of the shape measurement.
There are two main reasons for performing the flat calibration: (i) determination of the illumination angle change dθ for every measurement step that is essential in order to retrieve depth information, see Eq. (7), and (ii) measurement of the phase difference maps dP TLT used for background removal [43].
As aforementioned, two source point digital holographic contouring exhibits a hyperboloid background; however, it can be approximated by a tilted plane. Such background (corresponding to dP TLT in Eq. (6)) can be measured during the flat calibration and subtracted when a sample under test is measured. Moreover, such background removal enables processing of the data without the application of any spatial phase unwrapping algorithm since the high spatial frequency due to the illumination angle variation is removed.
Nine digital holograms were captured at different M5 mirror positions (different illumination angles) as it rotated. The measured phase fields at the nine different mirror positions were processed using Eq. (11). This results in a sequence of phase differences d P of which some are shown in Fig. 3(a). The numbers in the upper right corners denote for what mirror positions the phase difference map is shown (e.g. Mirror position 4-3 denotes the phase difference between the 4 th and 3 rd measurements). The sequence of the phase difference arrays in Fig. 3(a) was saved and later used for the background removal. As clear from Fig. 3(a), even for the temporal (differential) processing using two subsequent positions, spatial unwrapping is needed. However, unlike the TM (or other samples to be measured) the flat sample can be made of any material having optimal roughness and reflectivity so the high phase quality can be reached for ease of unwrapping. Increasing the number of measurement positions during the rotary stage movement could avoid spatial unwrapping for the flat sample. However, the Goldstein algorithm [44] was applied to spatially unwrap the phase maps sequence in Fig. 3(a) in this particular case. Further, the phase difference sequence d P was processed using Eq. (13) resulting in a sequence of unwrapped and cumulatively summed phase difference maps dP. The phase difference maps dP describe how the phase has changed for a mirror position in the upper right corner with respect to the first (initial) position of the mirror, see Fig. 3(b).
Since the sample is flat, for which z(x) = 0, Eq. (5) simplifies to: so the measured phase dP represents the background dP TLT . The illumination beam angle change dθ can be determined from Eq. (15) as a slope of the phase map with respect to the x-axis: In Eq. (16), the slope ∂dP⁄∂x was obtained as a parameter of 1st order polynomial fitting by means of least square: dP=∂dP/∂x x+∂dP/∂y y + C, where C is an intercept. The polynomial fitting was applied on the measured phase maps dP shown in Fig. 3(b) and results of the fitting are plotted in Fig. 4. The fitting is reliable since dispersion of the coefficient ∂dP⁄∂x estimated as 99% confidence bounds ranges from 4×10 −6 to 10×10 −6 that is negligible when comparing to values of ∂dP⁄∂x ranging from 3×10 −3 to 74×10 −3 . The measured value of the maximal angle change 0.072°is very close to the estimation 0.075°based on the yardstick measurement explained above. Illumination angle change dθ as a function of mirror position is plotted in Fig. 5. There are two important points to note. First, the flat calibration is not carried out simultaneously with the sample measurement and therefore repeatability of the mirror M5 motion plays an important role. The introduced flat calibration is repeated five times when the flat calibration sample is removed and replaced back between the measurements. Error bars in Fig. 5 plot represent standard deviations of the repeated measurements. The average value of the standard deviations for all mirrors positions is σ(dθ) = 0.001°. Three times of the standard deviation value is used in section 3.2 to estimate the measurement error budget.
Second, we investigated the background phase map variation dP when the flat sample is placed at different z positions within the same objective (OBJ) depth of field. Standard deviation between dP maps at different z positions was within 0.08 radians. Using Eq. (9), such variation of dP introduces uncertainty of the depth of the object (z) within 10 µm. The depth of the TM is approximately 1.5 mm and therefore the uncertainty 10 µm introduces relative error about 0.7%. The flat calibration is therefore valid and robust.

Step-height sample measurements
After the flat sample calibration, background phase difference maps dP TLT (dθ) and values of illumination angle change dθ for every mirror position were saved. Before proceeding to TM shape measurement, the performance of the measurement apparatus was further verified on a stepped gauge provided by National Institute of Standards and Technology (NIST). The NIST sample has precisely defined height steps and its geometry and the measured area are illustrated in Fig. 6. The intensity data were captured and phase difference maps were retrieved. The phase sequence after the background removal dP is shown in Fig. 7(a). If the background phase difference maps dP TLT (dθ) and the sample phase maps have different tilts due to an improper placement (significant tilt) of the calibration flat sample, the phase difference maps dP TLT (dθ) can be numerically adjusted by adding/subtracting the numerically generated tilted plane [45] (not needed in our case). The shape z was computed as a slope dP/dk independently in a pixel-wise manner. The phase difference dP as a function of sensitivity factor dk at three representative points (A,B,C) corresponding to different height steps is plotted in Fig. 7(b) while the resulting shape is shown in Fig. 7(c). The surface shape as a wireframe mesh and 1D slice is displayed in Fig. 8(a) and Fig. 8(b).
In order to estimate the accuracy of the measurement, the black dashed area in Fig. 7(c) was used for further analysis. The histogram in Fig. 9(a) introduces representation of the distribution of measured shape data. The values are normally distributed about three significant peaks corresponding to individual steps of the NIST gauge, see Fig. 6. The NIST gauge surface is curved along x dimension but the maximal surface depth within the dashed area is approximately 10 microns. This value is much lower than expected uncertainty and thus does not influence significantly the results.   A Gaussian function representing the normal distribution was fit to the three parts of the histogram in order to estimate some statistical features, see red line in Fig. 9(a). In Eq. (17) the symbol σ stands for standard deviation and z 0 is the mean of the distribution. About 68% of the measured values lie within one standard deviation, 95% of the measured values within two standard deviations and 99.7% are within three standard deviations. The standard deviation σ, therefore, can be used for accuracy estimation, see Table 1.
The maximal angle span of dθ determines the accuracy in a way that the larger angle change the better sensitivity (accuracy) of the measurement. In Fig. 9(b), the standard deviation σ (within areas denoted by B) is plotted as a function of angle change dθ exhibiting a decrease in measurement uncertainty with increasing dθ.

Tympanic membrane measurement
A fresh non-fixed human postmortem temporal bone was used for the measurement. The cartilaginous and bony parts of the ear canal were removed to provide optical access to the TM. Unlike the multiple wavelength technique employing a low power tunable laser, the multiple angle shape measurement technique does not require any TM surface manipulations like spraying white paint to increase its reflectivity. The human temporal bone study was approved by the Institutional Review Board of the Massachusetts Eye and Ear Infirmary. The use of de-identified human cadaveric tissue (only age and gender information are available) does not require a human subject study protocol. The laser light intensity irradiating the TM was 127 W/m 2 while the exposure time was about 1s. For such low intensities the tissue is not thermally influenced and there are supposed to be no damages [46]. The maps in Fig. 10(a) represent phase map differences dP at different mirror positions mathematically described in Eq. (5) with respect to the first mirror position. As the illumination angle increases (higher number for mirror position), the number of interference fringes proportionally increases and for large angle change, spatial unwrapping can no longer be reliable due to low contrast of the dense fringes. On the other hand, the introduced differential approach using two consecutive measurements -mirror positions (see Fig. 10(b)) and background removal avoids spatial unwrapping. Considering the angle change between two consecutive phase difference maps dθ≈0.01°, the surface contours (see Eq. (10)) are distant ∆z=5.8 mm. The surface contours distance is larger than the maximal height of TM (typically about 1.5 mm) and thus no spatial unwrapping is needed. The sequence of the phase difference maps was processed by Eq. (13) and the shape (Fig. 11(a)) was retrieved as a slope dP/dk, see Eq. (14). The resolved TM shape in Fig. 11(a) can be used to understanding mechanics of the TM or diagnose middle ear pathologies. Moreover, the TM shape will be used for correction of displacement measurements [36]. The correction required surface normal vectors along the TM surface. However, the computation of surface normal vectors is sensitive to high spatial frequency errors. In order to separate the TM shape from the high spatial frequencies generated mainly by speckle character due to the diffusive surface measurement, an 11th order polynomial surface was fitted to the measured shape. The fitted surface shape is shown in Fig. 11(b) and Fig. 12 as 2D plot, 3D plot, respectively. A histogram of the residuals of the fitting (difference between the shapes in Fig. 11(a) and Fig. 11(b)) is plotted in Fig. 11(c). It exhibits the normal distribution character with the standard deviation σ=0.141 mm that gives a clue about the noise level. In addition, the 1D plot in Fig. 12(b) represents values before (red line) and after (blue line) the polynomial fitting. There are several contributors to the shape measurement error budget ∆z. The uncertainty of the wavelength ∆λ is lower than other sources like phase noise ∆dP, angle uncertainty ∆dθ, and uncertainty ∆θ due to yardstick measurement for the determination of the illumination angle θ. Therefore, the uncertainty ∆λ is not considered. Although the error sources can be mutually compensated, for the error estimation we can apply total differential: The phase error, experimentally determined as ∆dP=0.8 rad, is a sum of different distortion sources e.g. electronic noise, speckle decorrelation, background intensity variation, speckle noise or intensity ratio of reference and object wave, diffraction patterns from dust particles, environmental distortions, etc. The value was estimated from a standard deviation map obtained from the consecutive phase fields, see for instance Fig. 10(b). The median of three times the standard deviation map was used as the uncertainty value ∆dP . The uncertainty due to the angle change ∆dθ=0.003°can be estimated as three standard deviations of the repeatability of the mirror positioning graphically introduced as error bars in Fig. 5. The illumination angle θ was computed using a triangulation considering length measurement by a tape measure resulting in the uncertainty estimation ∆θ=1.5°. Error contribution from the individual sources (computed as terms in Eq. (18)) as well as their root-mean-square (total error) are summarized in Table 2.

Future work -concept of a digital holographic otoscope
Within the research team, a parallel development is ongoing using fiber optics for illumination as well as for imaging. The combination of the methodology introduced in this paper with fiber optics can result in a compact instrument (digital holographic otoscope) enabling in-vivo measurement of the TM through the intact ear canal. Apart from visual inspection performed by Otologists to assess the TM health, the device in development has a potential to measure steady (shape) as well as dynamic (deformation) properties of TM with only one instrument. The principal arrangement of the developing instrument is shown in Fig. 13. The average human ear canal is about 25 mm in length and 7-9 mm in diameter. A working distance of 10 -15 mm between the instrument face (speculum-like shaped) and the TM is considered. The imaging system will have an angular field of view of approximately 60 degrees so that the whole TM surface can be imaged. The front lens is planned to be about 2 mm in diameter, and will allow space for illumination optics within the speculum. The illumination of the sample will be realized using a bundle of base fibers of 100 µm diameter. Either a fiber stretchers or acusto-optical modulators in heterodyne mode [47] will provide sufficiently rapid phase shifting. The sweeping of the angle dθ will be realized by switching the light between individual fibers using an electronically controlled optical switch. Therefore, no movable part will be required and the measurement can be performed quickly. Assuming 10 illumination directions, each taking 1ms, the measurement time cost is tens of ms (much more rapid than slow motions of the subject and TM) and thus makes quantitative measurements on live subjects possible. In order to increase illumination efficiency in terms of both wavefront and intensity distribution, all the fiber ends will be functionalized by printed optical elements. The worst case peak intensity impinging the TM is planned to be 150 W/m 2 for less than 1 second exposure time that safely meets limiting value of irradiation dose of tissue according to IEC 60825. The illumination part and the imaging part will be placed on opposite sides of the instrument face, see Fig. 13 (front view). Due to limited diameter of the ear canal, the illumination angle θ is limited to approximately 15°. As follows from Eq. (7), the smaller illumination angle results in lower sensitivity of the method. The sensitivity is roughly half when compared to the experiments on the optical bench introduced in section 3.2 where θ∼30°. On the other hand, assuming a distance span between the most apart bundled fibers 1 mm, the sweeping angle for the working distance will be more than 3°and thus the lower sensitivity due to smaller illumination angle θ will be compensated by the larger sweeping angle dθ. Therefore, we expect the error budget not to be worse than values introduced in Table 2. The data processing will be very similar to the evaluation described above, however, the calibration procedure may be slightly different.

Conclusions
This study provides a new method for the shape measurement of the tympanic membrane. The method is based on multiple angle digital holography employing a single wavelength laser. Using this principle, it is possible to compute the complete surface topology map of the TM without the need for spatial unwrapping of phase difference maps.
We proposed and constructed the experimental setup, which was used to demonstrate the measurement principle and its functionality. The method was successfully verified by measuring a reference stepped gauge provided by NIST. The uncertainty of the method was estimated to be several tens of microns and there is a room for improvement by a simple change of the arrangement's geometry. In order to demonstrate the ultimate possibilities of the developed method, a fresh non-fixed human postmortem TM was used as a test object. Potential application of the proposed methodology for in vivo measurements in live human ears is discussed, which may lead to development of an objective diagnostic tool to help differentiate middle-ear diseases in otology and audiology clinics.
Generally, the method is not limited to TM shape measurements but is also suitable for characterization of objects with complex geometries made of common materials (such as metals, plastics, etc.), as well as for the characterization of complex composite structures such as acoustic metamaterials, active acoustic metasurfaces, etc. The method can also be adopted into various lateral scales ranging from digital holographic microscopy using microscopic objectives to the shape measurements of large and distant surfaces.