ABSTRACT Title of Dissertation: DEVELOPMENT AND APPLICATION OF COMPUTER-AIDED FRINGE ANALYSIS Zhaoyang Wang, Doctor of Philosophy, 2003 Dissertation directed by: Associate Professor Bongtae Han Department of Mechanical Engineering Photo-mechanics methods have matured and emerged as important engineering tools. Although numerous image-processing algorithms have been developed to complement interferometric measurement techniques, these algorithms have been implemented originally for classical interferometry and the extensions to the general photo-mechanics fringe analysis have been limited. In this dissertation, the existing computer-aided digital fringe image analysis and processing techniques are investigated; the most appropriate fringe image processing schemes and their limitations are identified. To make the computer-aided fringe analysis practical to the real engineering problems, the existing schemes are improved and a series of new fringe analysis techniques are developed. Among these new techniques, the self-adaptive fringe filtering scheme considers not only the orientations of the local fringes but also the local fringe densities; the enhanced random phase shifting algorithm can detect the phase shift amounts and the full-field phase distributions automatically and simultaneously; the hybrid semi- automatic O/DFM fringe centering technique combines the advantages of existing techniques and can be employed to obtain full-field fractional fringe orders and their gradients accurately. Based on the study, a Windows GUI-based expert software system is developed for interferogram fringe analysis and processing. This expert system includes all the algorithms presented in this dissertation. Selected but original applications of the computer-aided fringe analysis are presented. They include: (1) development of infrared diffraction interferometer for co- planarity of high-density solder bump patterns; the infrared light enables the regularly spaced solder bump arrays to produce well-defined diffracted wavefronts, (2) development of an inverse method to determine elastic constants using circular disc and moir? interferometry; this method uses a non-linear over-deterministic approach to determine elastic constants simultaneously, and (3) applications to out-of-plane shape and warpage measurement and in-plane displacement and strain measurements of electronic packaging components. DEVELOPMENT AND APPLICATION OF COMPUTER-AIDED FRINGE ANALYSIS by Zhaoyang Wang Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2003 Advisory Committee: Associate Professor Bongtae Han, Chair Professor Donald Barker Assistant Professor Hugh Bruck Professor Jaime C?rdenas-Garc?a Professor Rama Chellappa ? Copyright by Zhaoyang Wang 2003 DEDICATION To my parents and my wife for their support and encouragement in this endeavor ACKNOWLEDGEMENTS I would like to express my sincere gratitude and appreciation to my advisor, Prof. Bongtae Han, for his excellent guidance, support and enthusiastic encouragement throughout this work. I am indebted to Prof. Jaime Cardenas-Garcia for his generous assistance in the circular disc compression experiment. I am also very grateful to Prof. Hugh Bruck for the useful discussions, to Prof. Donald Barker and Prof. Rama Chellappa for serving on the advisory committee. Finally, thanks are also extended to Mr. Michael Mello at Intel corporation and my friends and colleagues here at the University of Maryland, especially Sungyeol Cho and Seungmin Cho, for their assistance. i TABLE OF CONTENTS List of tables .................................................................................................v List of figures...............................................................................................vi 1. Introduction and background................................................................. 1 1.1 Problem statement............................................................................................1 1.2 A review of digital fringe analysis and processing............................................7 1.2.1 Intensity-based analysis................................................................................7 1.2.2 Phase measurement technique ......................................................................8 1.3 Scope and objective of this dissertation............................................................9 2. Pre-processing of fringe images ...........................................................11 2.1 Spatial averaging and low-pass filtering .........................................................11 2.2 Spatial median filtering ..................................................................................13 2.3 Frequency filtering: Fourier transform low-pass filtering ...............................14 2.4 Limitation of the conventional filtering methods ............................................15 2.5 New self-adaptive filtering.............................................................................15 2.5.1 Orientational filtering.................................................................................16 2.5.2 Selection of filtering window based on fringe density.................................17 2.5.3 New self-adaptive filtering .........................................................................18 ii 3. Fringe analysis I: Automatic analysis ...................................................22 3.1 Single image: Fourier transform method .......................................................22 3.1.1 Principle of general Fourier transform method ...........................................22 3.1.2 Principle of carrier Fourier transform method.............................................25 3.1.3 Discussions on Fourier transform method...................................................29 3.2 Multiple images: Phase shifting method.........................................................33 3.2.1 Principle of phase shifting method .............................................................33 3.2.2 Enhanced random phase shifting algorithm ................................................38 3.3 Fractional fringe order calculation: phase unwrapping....................................44 3.3.1 Sequential filling........................................................................................45 3.3.2 Minimum spanning tree..............................................................................46 3.3.3 Preconditioned-conjugate-gradient (PCG) least-square iteration.................46 3.3.4 About the phase unwrapping methods ........................................................48 3.4 Calculation of fringe order gradient................................................................50 3.5 Limitations of automatic analysis ...................................................................52 4. Fringe analysis II: Semi-automatic analysis..........................................53 4.1 Single image: fringe centering method ...........................................................53 4.1.1 Fringe centerline detection .........................................................................54 4.1.2 Semi-automatic fringe orders assignment ...................................................60 4.1.3 Fractional fringe order calculation: fringe order interpolation.....................61 4.1.4 Fringe order gradient or strain calculation ..................................................66 4.2 Overall limitations of fringe centering method ...............................................66 iii 4.3 New hybrid O/DFM fringe centering method for multiple images..................69 4.3.1 Optical/digital Fringe multiplication (O/DFM) ...........................................69 4.3.2 Hybrid approach.........................................................................................72 4.3.3 Advantages and disadvantages of hybrid semi-automatic analysis ..............73 5. Software development..........................................................................74 6. Application I: Out-of-plane shape and warpage measurement ..............81 6.1 Infrared diffraction interferometer for co-planarity measurement of high density solder bump pattern...................................................................................................81 6.1.1 Introduction ...............................................................................................82 6.1.2 Infrared diffraction Fizeau interferometry...................................................88 6.1.3 Co-planarity of flip-chip package bump pattern..........................................96 6.1.4 Summary..................................................................................................101 6.2 Warpage measurement of electronic packaging components ........................101 6.2.1 Warpage measurement using Twyman-Green interferometry ...................102 6.2.2 Warpage measurement using shadow moir? method.................................105 6.2.3 Warpage measurement using far infrared Fizeau interferometry (FIFI).....108 7. Application II: In-plane displacement and strain measurement...........111 7.1 Inverse method to determine elastic constants using circular disc and moir? interferometry..........................................................................................................111 7.1.1 Introduction .............................................................................................111 7.1.2 Background..............................................................................................113 iv 7.1.3 Over-deterministic inverse approach ........................................................118 7.1.4 Technical considerations for implementation............................................123 7.1.5 Experimental verification.........................................................................131 7.1.6 Discussion................................................................................................137 7.1.7 Summary..................................................................................................138 7.2 In-plane deformation and strain measurement of electronic packaging components .............................................................................................................139 7.2.1 In-plane deformation measurement using moir? interferometry................139 7.2.2 In-plane deformation and strain measurement using microscopic moir? interferometry......................................................................................................142 7.2.3 Coefficient of thermal expansion (CTE) measurement using moir? interferometry......................................................................................................147 8. Conclusions........................................................................................150 8.1 Conclusions .................................................................................................150 8.2 Future work .................................................................................................153 References ................................................................................................155 v LIST OF TABLES Table 1.1 In-plane photomechanics methods ...................................................................3 Table 1.2 Out-of-plane photomechanics methods ............................................................4 Table 7.1 Experimental results ....................................................................................137 Table 7.2 Results of CTE measurement at 60~100 ?C..................................................149 vi LIST OF FIGURES Figure 1.1 Illustration of large random noise. ..................................................................6 Figure 2.1 Example of low-pass filtering .......................................................................12 Figure 2.2 Example of median filtering .........................................................................13 Figure 2.3 Example of Fourier transform low-pass filtering ...........................................14 Figure 2.4 Schematic diagrams of the local fringe direction...........................................16 Figure 2.5 Example of fringe image filtering .................................................................21 Figure 2.6 Example of self-adaptive filtering window....................................................21 Figure 3.1 Example of fringe analysis using general Fourier transform technique..........24 Figure 3.2 3-D mapping using Fourier transform technique ...........................................28 Figure 3.3 Processing ideal fringes with in-side edges using Fourier transform..............31 Figure 3.4 Processing ideal fringes with incorrect shift of first harmonic .......................33 Figure 3.5 Example of 5-frame phase-shifting algorithm ...............................................36 Figure 3.6 Peak-to-valley (P-V) phase error versus percent phase shifter error [61] .......39 Figure 3.7 Example of a random phase shifting .............................................................43 Figure 3.8 Wrapped and unwrapped phase for a linearly changing displacement field ...44 Figure 3.9 Unwrapping of complicated phase using different algorithms .......................49 Figure 3.10 Examples of strain calculation using different gage lengths.........................51 Figure 4.1 Flow chart of fringe centering method ..........................................................53 vii Figure 4.2 Pixel matrix of 5?5 and directions for fringe peak detection .........................55 Figure 4.3 Example of fringe centerline detection..........................................................59 Figure 4.4 Example of semi-automatically fringe order assigning..................................60 Figure 4.5 1-D interpolation ..........................................................................................62 Figure 4.6 Improved 1-D interpolation ..........................................................................64 Figure 4.7 Example of improved 1-D interpolation........................................................64 Figure 4.8 Example of 2-D interpolation........................................................................65 Figure 4.9 Example of strain calculation........................................................................66 Figure 4.10 Example of fringe analysis using fringe centering method ..........................68 Figure 4.11 Schematic illustration of O/DFM ................................................................70 Figure 4.12 Example of O/DFM....................................................................................71 Figure 4.13 Example of strain calculation......................................................................72 Figure 5.1 Schematic architecture of the fringe analysis methods in software ................78 Figure 5.2 Examples of screen captures of the software.................................................80 Figure 6.1 Optical configuration of classical interferometry...........................................86 Figure 6.2 Schematic illustration of reflection from the bump arrays and the underlying spaces .........................................................................................................87 Figure 6.3 Basic concept of IR diffraction interferometry..............................................88 Figure 6.4 Required tilt angle for different bump pitches when the first diffraction order is used.........................................................................................................90 Figure 6.5 Changes of optical path length when point P moves to P?..............................91 Figure 6.6 Optical configuration of IR diffraction interferometry...................................94 viii Figure 6.7 Fringe patterns obtained from (a) the conventional IR Fizeau interferometry and (b) the proposed method.......................................................................96 Figure 6.8 Four phase-shifted images obtained from a flip-chip bump arrays on an organic substrate.........................................................................................97 Figure 6.9 Results of fringe processing..........................................................................98 Figure 6.10 Three dimensional representation of bump co-planarity..............................99 Figure 6.11 Deviations along (b) AA? and (c) BB? .......................................................100 Figure 6.12 Warpage measurement of a TAB package.................................................104 Figure 6.13 Principle of shadow moir? method............................................................106 Figure 6.14 Surface measurement of a quarter coin......................................................107 Figure 6.15 Warpage measurement of a PBGA package ..............................................108 Figure 6.16 Warpage measurement of a FC-PBGA package ........................................110 Figure 7.1 Schematic diagram of moir? interferometry................................................114 Figure 7.2 (a) Schematic diagram of a circular disc in compression, and (b) U and (c) V displacement fields obtained from a theoretical solution............................117 Figure 7.3 Sensitivity of displacement fields to elastic constants, E and ?....................119 Figure 7.4 (a) Theoretical fringe patterns with systematic errors caused by rigid-body displacements, and (b) the corresponding fringe patterns after adding random noise with a maximum variance of 0.25 fringe. .........................................124 Figure 7.5 Schematic diagram of the areas used in the nonlinear least squares calculation. .................................................................................................................126 Figure 7.6 Effect of calculation area on ? and ?. .........................................................127 ix Figure 7.7 Effect of location of the origin on ? and ?. .................................................129 Figure 7.8 (a) Schematic diagram of the specimen; (b) Schematic diagram of compression loading fixture used in the moir? experiments; (c) Photo of the experiment setup and (d) the specimen loaded in the loading fixture. ........133 Figure 7.9 U field fringe patterns obtained at 1800 N; (a)-(d) four phase-shifted moir? patterns; (e) and (f) phase maps before and after unwrapping process........134 Figure 7.10 U field fringe patterns obtained at 2200 N; (a)-(d) four phase-shifted moir? patterns; (e) and (f) phase maps before and after unwrapping process........135 Figure 7.11 Values of ? and ? obtained from the fringe patterns in Fig. 7.8 and 7.9.....136 Figure 7.12 Moir? fringe pattern of PQFP package with hygroscopic and thermal deformations.............................................................................................140 Figure 7.13 2-D deformations......................................................................................141 Figure 7.14 Schematic diagram of the flip-chip assembly on a high-density substrate (a) before and (b) after specimen preparation .................................................143 Figure 7.15 Micrographs of the region of interest. .......................................................144 Figure 7.16 Microscopic U and V displacement fields of (a) the bare substrate and (b) the flip-chip assembly.....................................................................................145 Figure 7.17 Deformations of the solder and the metal via in substrate..........................146 Figure 7.18 Deformations of the solder and the metal via in Flip-chip assembly. .........146 Figure 7.19 Strain ?y distributions in the solder and the metal via. ...............................147 Figure 7.20 CTE measurement of PCB........................................................................148 Figure 7.21 CTE measurement of module....................................................................149 1 1. INTRODUCTION AND BACKGROUND 1.1 Problem statement In recent years, various photomechanics methods have matured and emerged as important engineering tools. The methods provide the full field information of displacement field. Representative examples of the methods used for a quantitative deformation analysis are shown in Table 1.1 and Table 1.2 for in-plane and out-of-plane displacements, respectively [1]~[6]. The fringe patterns represent the contours of equal displacements. The intensity distribution of the fringe patterns can be expressed as )]y,x(cos[)y,x(I)y,x(I)y,x(I am ?+= (1.1) where I is the intensity distribution of the interferogram, Im is the mean intensity, Ia is the intensity modulation amplitude, ? is the angular phase information of the interferogram, and (x,y) represents all the points in the x-y plane of the object and the interferogram; ? represents the fringe order N at each point of the pattern by (x,y)2N(x,y)?=pi . Experimental fringe analysis is a procedure to obtain the desired experiment results or physical quantities from the fringe patterns. In general, only the integer fringe orders can be directly determined from the fringe patterns, which are located along the centerlines of the black fringes. The fringe analysis involves two steps: one is how to obtain the 2 fractional fringe orders and the other is how to get the desired results from the fractional fringe orders. The fringe patterns have been analyzed manually for experiment analyses. As digital image processing has become more accessible, automatic analyses of the fringe patterns have become practical [7]. Numerous image-processing algorithms have been developed to complement interferometric measurement techniques. The methods utilize a single or a series of phase-shifted interferograms to compute the fractional fringe orders. The algorithms were originally implemented for classical interferometry. Their extensions to general fringe analyses have been limited because of inherent optical noise encountered in the modern photomechanics techniques. Under an idealized condition where the interferogram is described faithfully by Equation (1.1), mathematical determination of the phase information can be readily achieved. If random noise with relatively small amplitude exists, the general image processing algorithms for averaging and smoothing can be utilized effectively to eliminate the noise. In practice, however, the general image processing algorithms are not directly applicable because of geometrical discontinuities in the region of interest as well as random noise with large amplitude. 3 Table 1.1 In-plane photomechanics methods Method Basic principle Sensitivity (fringes/disp.) Contour Interval (disp./fringe) Field of View Example Geometric moir? Additive intensities Less than 100 lines/mm Greater than 10 micrometers Large Moir? interferometry Laser interference 2.4 lines/micrometer 0.417 micrometers Small (typically 5mm to 50mm) 4 Method Basic principle Sensitivity (fringes/disp.) Contour Interval (disp./fringe) Field of View Example Microscopic moir? interferometry Interference under refractive- index medium and microscope 4.8 lines/micrometer 208 to 20.8 nanometers Microscopic (typically 50 micrometers to 1 mm) Table 1.2 Out-of-plane photomechanics methods Method Basic principle Sensitivity (Fringes/disp.) Contour Interval (disp./fringe) Field of View Example Shadow moir? Additive intensities Less than 100 lines/mm Greater than 10 micrometers Large (up to 100 mm) 5 Method Basic principle Sensitivity (Fringes/disp.) Contour Interval (disp./fringe) Field of View Example Infrared Fizeau interferometry Light interference 200 to 400 lines/mm 2.5 to 5 micrometers Medium (5 to 45 mm) Twyman- Green interferometry Laser-light interference 3.2 lines/micrometer .317 micrometers Small (less than 5 mm) 6 The inherent optical noise of an interferogram is illustrated in Figure 1.1. The fringe pattern represents a horizontal displacement of a layered metal/ceramic composite subjected to a uniform bending, which was obtained from a microscopic moir? experiment. The intensity distributions along the cross lines in the fringe patterns are also shown in the figure. In spite of well-defined black fringes with excellent contrast, random noise with extremely large amplitude is evident. The noise was caused by imperfection of diffraction gratings used in the experiment. In addition, a strong gradient (or strain) is expected at the interfaces of material discontinuity. If a general image- processing algorithm is employed, the original deformation field can be altered significantly. Image processing routines with specific functions are required to process the fringe pattern. Figure 1.1 Illustration of large random noise. 7 1.2 A review of digital fringe analysis and processing Digital fringe image processing techniques fall into two categories: one is an intensity-based analysis and the other one is a phase measurement technique. The former is a semi-automatic processing technique while the latter can be regarded as an automatic processing technique. The intensity-based fringe analysis was established based on the traditional manual processing method. This technique involves detecting and thinning fringe centerlines (skeletons), assigning fringe orders and interpolating fringe orders. The intensity-based fringe processing is a tedious procedure and development of this technique has not been active since 1990?s [8]~ [38] because of its inherent limitation. The phase measurement technique [39]~[97] offers an automatic processing procedure and has been widely used since late 1980?s. In spite of its popularity, the digital fringe image processing is still not mature for practical applications. One of the reasons is that the techniques were established based on theoretical descriptions of the optical interferograms; significant improvement is required to make the techniques practical for real engineering applications. 1.2.1 Intensity-based analysis The core of the intensity-based analysis technique [8]~[38] includes fringe centerline detection and fringe order interpolation. There are two kinds of methods for extracting the fringe centerlines from fringe patterns. One method involves binarizing the fringe 8 patterns and then skeletonizing the binary fringes; another method finds the fringe centerlines through detecting the local maxima and minima of fringe intensities. Compared with binarization method, the fringe peak detection method usually yields centerlines with less error; however, this method is more sensitive to noise. Although many fringe centerline detection methods have been proposed and developed, few study work can be found in the literature for automatic fringe order assignment and interpolation. Generally, it is practically impossible to assign fringe orders automatically because the fringe patterns do not contain the fringe order information. A typical solution to this problem is detecting the ascending or descending directions of the fringe order while changing the load in the experiment. For fringe order interpolation, only one-dimensional (1-D) algorithms based on linear interpolation have been implemented. It is obvious that 1-D interpolation is not adequate for most cases simply because the fringe patterns contains two-dimensional (2-D) information. 1.2.2 Phase measurement technique Unlike the intensity-based fringe analysis which utilizes only a small portion of the fringe pattern, i.e., the integer fringe orders, the phase measurement method uses full- field fringe information for the analysis. The phase measurement technique includes Fourier transform analysis and phase shifting analysis. Fourier transform method [39]~[56] requires a single fringe image for the analysis. In order to separate the pure phase information in the frequency domain, Fourier transform usually requires carrier fringes; this brings difficulty in practice 9 because the frequency of the carrier fringe must be controlled accurately. Another critical limitation of Fourier transform technique is inability of handling discontinuities. Phase shifting method [57]~[77] is a widely used automatic fringe processing technique. This method uses a series of phase-shifted interferograms to calculate the fractional fringe orders. Theoretically, the phase shifting algorithm is ideal for a general fringe analysis. However, in real applications, the phase shifting technique also suffers from localized inaccuracy, especially when the gradients of the fringe orders are to be determined. Another important issue associated with phase shifting technique is phase unwrapping [78]~[97] process. Numerous phase unwrapping algorithms are available and new phase unwrapping algorithms are being proposed. It is worth noting that most of the algorithms aim to solve one or more specific problems. An investigation of the existing algorithms for automatic fringe analysis is in high demand. 1.3 Scope and objective of this dissertation The scope of this dissertation invovles investigation of the digital fringe processing techniques for computer-aided fringe analyses. The ultimate objective of this dissertation is to develop a general-purpose computer-aided fringe analysis tool and apply it to advanced engineering problems. The specific goals include: (1) Full-field mapping of fractional fringe order, which usually denotes a deformed configuration, and (2) Full-field mapping of the gradient of the fractional fringe order, which usually denotes the displacement gradient, i.e., strain. 10 To achieve the goals, the following tasks will be implemented and completed in this dissertation: (1) Survey and investigate the existing fringe image processing schemes and identify the most appropriate image processing schemes and their limitations; (2) Improve the existing schemes and develop new techniques to cope with the limitations of existing schemes; (3) Develop an expert software system for computer-aided fringe analysis and processing; (4) Apply the computer-aided fringe analysis expert system to the real engineering problems. 11 2. PRE-PROCESSING OF FRINGE IMAGES This chapter is devoted to the algorithms for pre-processing of fringe patterns to eliminate the noise. The existing algorithms are reviewed and the limitations are discussed. New hybrid algorithms are proposed for the fringe patterns and the results from the proposed algorithms are presented. 2.1 Spatial averaging and low-pass filtering Low pass filtering, also known as ?smoothing?, was developed to remove noise with high spatial frequency [7] ~[11]. This type of noise is often produced during the analog- to-digital conversion process (physical conversion of light energy into electrical signal). A typical form of low-pass filters is a moving window operator. The operator affects one pixel at a time, changing its value by some function of a ?local? region of pixels (?covered? by the window). The operator ?moves? in the x and y direction to change the entire image. The most common low-pass filter is a neighborhood-averaging filter. The neighborhood-averaging filters replace the value of each pixel, say ( )y,xI , by a weighted-average of the pixels in some neighborhood around it, i.e. a weighted sum of ( )qy,pxI ++ , with p = -k to k, q = -k to k, where k is a positive integer; the weights are given to each pixel in the window, usually non-negative value. If all the weights are equal, it is called a ?mean filter?. Mathematically the neighborhood-averaging filtering method can be expressed as: 12 ( ) ( ) ( )? ? ?= ?= ?++=? k kp k kq q,phqy,pxIy,xI (2.1) where ( )q,ph is the weight of averaging at pixel (p, q). Figure 2.1 shows a horizontal displacement fringe pattern of a layered metal/ceramic composite subjected to a uniform loading; Figure 2.1 (a) is the original fringe pattern and Figure 2.1 (b) is the same fringe pattern after low-pass filtering using a 5?5 window mean filter. Significant reduction of noise is evident. The neighborhood-averaging filter has an excellent smoothing effect. For the same reason, it can blur the images, especially along the edges, and thus it should be used with a caution when a large random noise is present. (a) Original fringe pattern (b) Fringe pattern after filtering Figure 2.1 Example of low-pass filtering 13 2.2 Spatial median filtering Median filtering [7]~[11] is a non-linear filtering method, and it is very useful for reduction of salt-and-pepper type noise (i.e. isolated noise with extreme values). The median filtering replaces each pixel value by the median of its neighbors. For instance, consider a 3 by 3 pixel window, the intensity values of the nine pixels in the window are sorted in an ascending order and the value of the pixel in the center of the window is replaced by the fifth largest value. Figure 2.2 shows the result of applying a 5?5 window median filtering on the same image in Figure 2.1 (a). Figure 2.2 Example of median filtering In general, median filtering takes longer time to execute because of a logical operation required in the algorithm and it is not very effective for noise reduction of fringe patterns. However, the method does not alter the edges as much as low pass 14 filtering. It is very effective for the filtering of phase map, which contains pixels with undefined phase values. The phase map issue will be discussed more in Chapter 3. 2.3 Frequency filtering: Fourier transform low-pass filtering Fourier transform low-pass filtering method [11] is a form of low pass filtering in frequency domain. With this method, an image is transformed into a frequency domain and an appropriate window is selected in the frequency domain to eliminate the high frequency noise. The filtering is expressed as: ( ) ( )[ ] ( ){ }v,uHy,xfFTFTy,xg 1 ?= ? (2.2) Figure 2.3 Example of Fourier transform low-pass filtering where FT represents Fourier transform, FT-1 represents inverse Fourier transform, H(u, v) is a filter function, f(x, y) is the original intensity and g(x, y) is the modified intensity. 15 The typical filter function includes Ideal filter, Butterworth filter, Exponential filter and Trapezoidal filter [11]. Figure 2.3 shows the result of the same fringe pattern in Figure 2.1 (a) after applying a Fourier transform Exponential filtering. 2.4 Limitation of the conventional filtering methods The spatial image filtering methods were developed for general photographic images and they are not most effective for fringe patterns, which have more distinct orientational intensity variations compared with the general photographic images. More specifically, the intensity of fringe patterns changes much more rapidly along the orientation perpendicular to the fringe lines than along the orientation parallel to the fringe. Furthermore, the spatial filtering methods do not consider the variation of fringe density. Consequently, a fixed size of conventional square window is not most desirable for fringe patterns. The similar limitations can be realized with Fourier transform filtering. The Fourier transform method filters the fringe in the global frequency domain. It considers neither the variation of fringe density nor the fringe orientation. 2.5 New self-adaptive filtering A new self-adaptive filtering method is proposed in this dissertation. The self- adaptive filtering is based on spatial averaging filtering; however, it detects both the orientations and the densities of the fringes during the filtering. 16 2.5.1 Orientational filtering Orientational filtering utilizes the essentially ?line? characteristic of the fringes. In the process, the fringe orientation at every pixel is determined first. Generally, a proper- size window (usually a square window) is used to calculate the variances of fringe intensities in each orientation of the window: ? ? = = ? ? ?? ? ? ?= k 1i 2k 1i j i j i j I k 1IVar (2.3) where I denotes the fringe intensity or pixel gray level, the superscript ?j ? denotes the jth fringe orientation, the subscript ?i? denotes the ith pixel in each orientation and k is the number of pixels in each orientation. Figure 2.4 Schematic diagrams of the local fringe direction 1 2 3 4 5 6 7 8 Orientation 1 Orientation 2 Fringe orientation definition Orientation 3 Orientation 4 17 Figure 2.4 illustrates an example. The example shows 8 possible orientations within the window. Details of the first four orientations are also shown for the window of 11?11. Among all the orientations, the orientation with the minimum fringe intensity variance is chosen as the fringe orientation. Then, instead of a square window which does not consider the fringe orientation, a rectangular window mask in the orientation of the fringe is used to smooth the fringes. 2.5.2 Selection of filtering window based on fringe density In the neighborhood-averaging filtering, the filters can be applied repeatedly or a large window size can be used to tailor the degree of noise smoothing. However, excessive smoothing can reduce the contrast of high-density fringes significantly, and thus can disturb the original phase information of the fringes. To cope with the problem, the filtering window size can be selected automatically based on the local fringe density. To find the local fringe density, firstly, the original fringe image is pre-smoothed using a small window neighborhood-averaging filtering. Then, the fringe image is binarized to be a monochrome (black-white) image with a threshold which can be the mean intensity of the fringe pattern. From the monochrome fringe image, the local fringe density along horizontal direction fx and along vertical direction fy can be determined; and the total local fringe density is calculated by 2 y 2 x fff += (2.4) 18 Finally, the original fringe pattern can be smoothed again using a low-pass filter with varying window whose size matches to the local fringe density. 2.5.3 New self-adaptive filtering In the proposed self-adaptive filtering, the orientational filtering and the fringe- density-based window filtering are combined to provide an optimum spatial filtering condition. The filter window is a rectangular in the orientation of the local fringes; the width and height of the rectangular window can be selected to be one pitch of the local fringes (i.e., reciprocal of fringe density) and one quarter of the pitch, respectively. In the real application, the detection of fringe orientation at every pixel can be erroneous because of the noise. A solution to this problem is to segment the whole fringe image into numerous blocks so that the fringe density and orientation are reasonably uniform in each block. Figure 2.5 illustrates the effect of different filtering methods. The original image is shown in (a), which represents a vertical displacement field of a notched specimen subjected to a cyclic loading, obtained by moir? interferometry. Using a conventional averaging filter with 3x3 size window, the noise in the high-density fringes around the crack was eliminated effectively but the noise in the low-density fringes at a far field was not cancelled (b). When a larger (7x7) window was used, the noise of the far field reduced significantly but it caused smearing of the high-density fringes (c). The result obtained from the proposed self-adaptive method is shown in (d) and effective reduction of noise in both areas is evident. The intensity distributions along the same white solid 19 line in (a) are shown in (e)~(h). Figure 2.6 illustrates the self-adaptive window at point A and B used in the filtering. The proposed algorithm eliminated the noise very effectively without altering the original shape of the intensity distribution. (a) Original image (b) 3?3 window averaging filtering (c) 7?7 window averaging filtering (d) Self-adaptive filtering 20 0 50 100 150 200 250 300 0 50 100 150 200 250 300 350 400 450 Location along the line (pixel) Pixel gray level (e) Intensity plot of (a) along the line 0 50 100 150 200 250 300 0 50 100 150 200 250 300 350 400 450 Location along the line (pixel) Pixel gray level (f) Intensity plot of (b) along the line 0 50 100 150 200 250 300 0 50 100 150 200 250 300 350 400 450 Location along the line (pixel) Pixel gray level (g) Intensity plot of (c) along the line 21 0 50 100 150 200 250 300 0 50 100 150 200 250 300 350 400 450 Location along the line (pixel) Pixel gray level (h) Intensity plot of (d) along the line Figure 2.5 Example of fringe image filtering 5x1, 45? window at A 45x13, 0? window at B Figure 2.6 Example of self-adaptive filtering window 22 3. FRINGE ANALYSIS I: AUTOMATIC ANALYSIS Two algorithms for automatic fringe analysis are studied: Fourier transform method [39]~[56] with a single image and phase shifting technique [57]~[77] with multiple images. The purpose of these algorithms is to find a phase information at every pixel, and thus to determine a fractional fringe order at every point in the fringe pattern. Mathematical descriptions are given and illustrations are followed. Limitations and restrictions in the application to the fringe patterns are discussed. 3.1 Single image: Fourier transform method 3.1.1 Principle of general Fourier transform method Repeated here, the intensity distribution of the fringe patterns can be expressed as )]y,x(cos[)y,x(I)y,x(I)y,x(I am ?+= (3.1) where I is the intensity distribution of the interferogram, Im is the mean intensity, Ia is the intensity modulation amplitude, ? is the angular phase information of the interferogram, and (x,y) represents all the points in the x-y plane of the object and the interferogram. The intensity function can also be expressed as ( ) ( ) ( ) ( )y,xcy,xcy,xIy,xI m ?++= (3.2) with ( ) ( ) ( )y,xia ey,xI21y,xc ?= (3.3) 23 where * denotes complex conjugate. After a two-dimensional discrete Fourier transform (DFT), the spatial frequency domain representation of the pattern becomes ( ) ( ) ( ) ( )hzhzhzhz ,C,C,A,F *++= (3.4) where A(?, ?) is the transform of Im(x, y), and C(?, ?) and C*(?, ?) are the positive and negative frequency spectra of the modulated carrier fringes. ? and ? are the spatial frequencies that represent intensity changes with respect to spatial distances. If the image size is 2n, the fast Fourier transform (FFT) can be used, which is much faster than the ordinary Fourier transform. At the frequency domain, if C(?, ?) can be isolated from A(?, ?) and C*(?, ?) in equation (3.4), then an inverse Fourier transform can be performed for C(?, ?). Finally, c(x, y) can be obtained at the spatial domain and the phase information ? can be calculated from ( )[ ] ( )[ ]y,xcRe y,xcImarctan=? (3.5) where Im[ ], Re[ ] represent the imaginary and real part of c(x, y), respectively. The phase ? obtained from the above equation ranges from -pi to +pi, and it does not reflect the fringe order. Phase unwrapping is required to make the phase represent the fractional fringe order and it will be discussed later in this chapter. In order to be able to isolate C(?, ?) from A(?, ?) and C*(?, ?) at the frequency domain, the intensity function should have continuous and monotonically changing derivatives across the field. Unfortunately, this condition is often violated for the fringe pattern representing engineering deformations. 24 (a) Original ideal fringes (b) Fourier spectra of (a) (c) Phase map from Fourier transform (d) Theoretical phase map Figure 3.1 Example of fringe analysis using general Fourier transform technique Figure 3.1 shows one example of fringe pattern analysis using general Fourier transform technique introduced above. The original computer-generated fringe pattern is shown in (a), The spectra of Fourier transform is shown in (b), where the positive first 25 harmonic (i.e., C(? ?) ) is isolated by a rectangular box shown in the picture. The result of the inverse Fourier transform of the isolated harmonic is shown in (c); it is a wrapped phase map. The corresponding theoretical phase map is shown in (d). The result illustrates that the above general Fourier transform method does not work for fringe analysis if the sign of the fringe gradient changes; i.e., it is difficult to isolate C(?, ?) at the frequency domain. A modified Fourier transform method was introduced to cope with this problem. 3.1.2 Principle of carrier Fourier transform method In this modified Fourier transform method, a spatial carrier fringe is added to the original pattern to make the sign of the fringe gradient unchanged across the fringe pattern. The carrier fringe is a uniform array of fringes representing a constant gradient. The original pattern is modulated by the relatively high frequency carrier pattern to produce a monotonically changing intensity function. This process is analogous to FM (frequency modulated) radio waves, where the information is the irregular part of an otherwise constant carrier frequency. The intensity function modulated by the carrier frequency can be expressed as ( ) ( ) ( ) ( ) ( )[ ]y,xyfxf2cosy,xIy,xIy,xI yxam ?++pi+= (3.6) where fx and fy are the linear components of the carrier frequency in the x and y direction, respectively. Again, the intensity function can be rewritten in the exponential form as ( ) ( ) ( ) ( ) ( ) ( )yfxfi2yfxfi2m yxyx ey,xcey,xcy,xIy,xI +pi??+pi ++= (3.7) with 26 ( ) ( ) ( )y,xia ey,xI21y,xc ?= (3.8) where * denotes complex conjugate. After a two-dimensional Fourier transform, the spatial frequency representation of the pattern becomes ( ) ( ) ( ) ( )yx*yx f,fCf,fC,A,F +?+?+????+??=?? (3.9) where A(?, ?) is the transform of Im(x, y), and C(?-fx, ?-fy) and C*(?+fx, ?+fy) are the positive and negative frequency spectra of the modulated carrier fringes. ? and ? are the spatial frequencies that represent intensity changes with respect to spatial distances. It is essential for the three terms in equation (3.9) to be completely isolated from one another in the frequency domain. This condition will be achieved if the signal of interest is sufficiently band-limited, and the frequency of the carrier pattern is large enough to separate the positive and negative spectra. Otherwise, the positive spectra and negative spectra will be overlapped in the frequency domain and the first harmonic cannot be separated as seen in the general Fourier transform method. At the frequency domain, C(?-fx, ?-fy) is isolated to eliminate A(?, ?) and C(?+fx, ?+fy) in equation (3.9). By shifting the center of the spectrum to the origin of the frequency axis, the carrier fx and fy are removed from C(?-fx, ?-fy). Then, the inverse Fourier transform is performed for C(?, ?), and finally c(x, y) can be obtained at the spatial domain. The phase information ? can be calculated from the same equation as equation (3.5), repeated here, ( )[ ] ( )[ ]y,xcRe y,xcImarctan=? (3.5) 27 Similarly, a phase unwrapping is required to make the phase represent the fractional fringe order. Because carrier Fourier transform can handle the practical engineering problems, it is more widely used than general Fourier transform method. The term Fourier transform actually represents carrier Fourier transform method. An example of the method is shown in Figure 3.2. The original fringe pattern is shown in (a), which contains two local maxima. The carrier frequency and the pattern modulated by the carrier are shown in (b) and (c), respectively. The modulated pattern has a monotonically changing gradient in the x direction. The result of the Fourier transform is shown in (d), where the positive first harmonic (i.e., C(?-fx, ?-fy) ) is isolated by a rectangular box shown in the picture. The first harmonic is shifted to the origin and the result of the inverse Fourier transform of the shifted harmonic is shown in (e); it is a wrapped phase map. The corresponding 3-D representation of the fractional fringe orders (or phase map) is shown in (f). (a) Original fringes (b) Carrier fringes 28 (c) Modulated fringes with carriers (d) Fourier spectra of (c) (e) Fringe phase map (f) 3-D unwrapped phase map Figure 3.2 3-D mapping using Fourier transform technique 29 3.1.3 Discussions on Fourier transform method The most important advantage of the Fourier transform method (or carrier Fourier transform method) is that the only one interferogram is required for the analysis. It is ideal for the automatic calculation of full-field displacements. Although ideal in the mathematical description, the Fourier transform method has several practical limitations. The most critical limitation is the lack of capability of handling discontinuities. As mentioned earlier, the DFT assumes that the original function has a cyclic variation with a continuous gradient. At the discontinuities, the process of transformation will spuriously distribute a large amount of energy over a wide range of frequencies. This not only obscures the signal power, but also makes isolation of the pure signal impractical. It is important to note that the edges of a fringe pattern are an inherent source of discontinuities even when the fringe pattern itself is continuous. The effect of discontinuities is illustrated in Figure 3.3, which represents an ideal U displacement field of a plate with a hole subjected to a uniform tension. The original and the modulated fringe patterns are shown in (a) and (b), respectively. The corresponding frequency spectra and the phase map obtained after the inverse transform are shown in (c) and (d). The theoretical phase map is shown in (e). As can be seen, the phase information along the edges is distorted significantly. 30 (a) Original fringe pattern (b) Modulated fringes with carriers (c) Fourier spectra of (b) (d) Phase map from Fourier transform 31 (e) Theoretical phase map Figure 3.3 Processing ideal fringes with in-side edges using Fourier transform It is also important to note that a precise shift of the first harmonic to the origin before the inverse transform is critical. The effect of precise shift is illustrated in Figure 3.4. (a) is the original ideal fringe pattern, the image size is 256 pixels by 256 pixels; (b) is the modulated fringes with carrier fx = 0.1 fringe/pixel, thus there are totally 25.6 fringes added through the whole horizontal image width; (c) is the corresponding spectra of (b). In this example, a precise shift of first harmonic should be 25.6 pixels; however, this is not allowed in digital image processing. (d) and (e) are the phase maps obtained from inverse Fourier transform after shifting the positive first harmonic to the left with 25 and 26 pixels, respectively; (f) is the correct phase map. The importance of the precise shift of the first harmonic is obvious in this example. In the real application, it is generally impossible to control the number of the carrier fringes to be exact integer number in the 32 digitized image. This is an inherent error from a digital processing and another practical limitation of the Fourier transform method. (a) Original fringe pattern (b) Modulated fringes with carriers (c) Fourier spectra of (b) (d) Phase with inadequate harmonic shift 33 (d) Phase with excessive harmonic shift (e) Correct phase map Figure 3.4 Processing ideal fringes with incorrect shift of first harmonic 3.2 Multiple images: Phase shifting method 3.2.1 Principle of phase shifting method The method utilizes a series of fringe or phase shifted interferograms to compute the fractional fringe orders [57]~[77]. The algorithms were originally implemented for classical interferometry. Recently, their applications have been extended for other advanced photomechanics methods. The basic algorithm and the enhanced algorithm are discussed. A new algorithm is proposed to cope with one of the most critical requirements of the method; namely, accurate phase shift. 34 3.2.1.1 3, 4, 5, N frames algorithm Repeated here, the intensity I(x, y) of an interferogram at a point (x, y) is given as [ ])y,x(cos)y,x(I)y,x(I)y,x(I am ?+= (3.1) There are three unknowns in the equation, namely Im, Ia, ?. Three simultaneous equations are needed to evaluate the unknowns. Experimentally, the three equations can be obtained by recording a series of intensity distributions with a uniform change of phase (or fringe order). The three equations can be expressed as ( ) ( ) ( ) ( )[ ]???+= y,xcosy,xIy,xIy,xI am1 (3.10) ( ) ( ) ( ) ( )[ ]y,xcosy,xIy,xIy,xI am2 ?+= (3.11) ( ) ( ) ( ) ( )[ ]?+?+= y,xcosy,xIy,xIy,xI am3 (3.12) where I1, I2 and I3 are the intensity distributions recorded with a phase change of ??, 0 and +?, respectively. From the three equations, the phase ?(x, y) can be determined as ( ) ( ) ( )( ) ( ) ( )? ? ? ?? ? ?? ?? ? ??=? y,xIy,xIy,xI2 y,xIy,xI sin cos1arctany,x 312 31 (3.13) When pi=? 32 , the above expression becomes ( ) ( ) ( )( ) ( ) ( )? ? ? ?? ? ?? ??=? y,xIy,xIy,xI2 y,xIy,xI3arctany,x 312 31 (3.14) For a more accurate phase calculation, other algorithms using more than three phase- shifted images have been developed. The most widely used algorithm uses four phase- shifted images. The set of four images are ( ) ( ) ( ) ( )[ ]y,xcosy,xIy,xIy,xI am1 ?+= (3.15) 35 ( ) ( ) ( ) ( ) ? ? ? ?? ? pi+?+= 2 1y,xcosy,xIy,xIy,xI am2 (3.16) ( ) ( ) ( ) ( )[ ]pi+?+= y,xcosy,xIy,xIy,xI am3 (3.17) ( ) ( ) ( ) ( ) ? ? ? ?? ? pi+?+= 2 3y,xcosy,xIy,xIy,xI am4 (3.18) Then, the phase is expressed in a simpler form as ( ) ( ) ( )( ) ( )? ? ? ?? ? ? ?=? y,xIy,xI y,xIy,xIarctany,x 31 24 (3.19) Another algorithm uses five images [61]. This algorithm was developed to minimize the cases of denominators with zero or near zero values, and thus to reduce uncertainties in the phase calculation. The five different intensities are obtained with a symmetric phase shift as ( ) ( ) ( ) ( )[ ]???+= 2y,xcosy,xIy,xIy,xI am1 (3.20) ( ) ( ) ( ) ( )[ ]???+= y,xcosy,xIy,xIy,xI am2 (3.21) ( ) ( ) ( ) ( )[ ]y,xcosy,xIy,xIy,xI am3 ?+= (3.22) ( ) ( ) ( ) ( )[ ]?+?+= y,xcosy,xIy,xIy,xI am4 (3.23) ( ) ( ) ( ) ( )[ ]?+?+= 2y,xcosy,xIy,xIy,xI am5 (3.24) The phase is given by ( ) ( ) ( )( ) ( ) ( )? ? ? ?? ? ?? ?? ? ??=? y,xIy,xIy,xI2 y,xIy,xI sin cos1arctany,x 513 42 (3.25) If pi=? 21 , 36 ( ) ( ) ( )( ) ( ) ( )? ? ? ?? ? ?? ?=? y,xIy,xIy,xI2 y,xIy,xIarctany,x 513 42 (3.26) Figure 3.5 Example of 5-frame phase-shifting algorithm With the phase shifting method, uncertainties in phase determination are present when the phase shifting amount, ?, is not correctly introduced during the measurement. 37 When the three-phase-step is used, the phase is more sensitive to the phase shift error. Small errors in the phase shift result in significant overestimation or underestimation of the phase. Using more phase steps, the errors in the phase shift can be smoothed out, which produces more stable results. Considering the computational time, four or five phase steps are most widely used in practice. Figure 3.5 shows an example of the 5- frame phase shifting algorithm. The last image is the calculated phase map. 3.2.1.2 Advanced phase shifting method: Carr? method With Carr? method [61], the phase shift amount is also treated as an unknown. The method uses four phase-shifted images as ( ) ( ) ( ) ( ) ? ? ? ?? ? ???+= 2 3y,xcosy,xIy,xIy,xI am1 (3.27) ( ) ( ) ( ) ( ) ? ? ? ?? ? ???+= 2 1y,xcosy,xIy,xIy,xI am2 (3.28) ( ) ( ) ( ) ( ) ? ? ? ?? ? ?+?+= 2 1y,xcosy,xIy,xIy,xI am3 (3.29) ( ) ( ) ( ) ( ) ? ? ? ?? ? ?+?+= 2 3y,xcosy,xIy,xIy,xI am4 (3.30) Assuming the phase shift is linear and does not change during the measurements, the amount of phase shift can be calculated as ( ) ( ) ( ) ( ) ??? ? ??? ? ?+? ???=? ? 4132 41321 IIII IIII3tan2 (3.31) and the phase at each point is determined as 38 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 14231 2314 142323141 2314 IIIItantan 2IIII IIII3IIIItan IIII ? ? ?????+???????= ?????? +?+?? ???? ?+????????????= +?+ (3.32) The advantage of Carr? algorithm is clear; it does not require accurate calibration of the phase shifting mechanism as long as it is linear and stable during the measurement. 3.2.2 Enhanced random phase shifting algorithm Under an idealized condition where the phase shifting is performed perfectly, any of the above algorithms would suffice to produce accurate phase information. In practice, the phase shifting has uncertainties. Figure 3.6 shows how each algorithm is sensitive to the phase shifting error, where peak to valley (P-V) phase errors are plotted as a function of the phase shift error. When a linear phase shift error (the actual phase shift error is a linear function of the desired phase shift) exists, the 3-frame and 4-frame algorithms yield a large error. The error is suppressed significantly with the 5-frame algorithm and it is nullified completely with the Carr? method. For a nonlinear phase-shifter error (the actual phase shift error is a non-linear function of the desired phase shift, and non-linear errors include non-linear response from a CCD detector, quantization error, vibration, etc.), however, all the methods yield substantial error in phase determination. 39 Figure 3.6 Peak-to-valley (P-V) phase error versus percent phase shifter error [61] It is important to note that the series of phase-shifted interferograms contain the information required to obtain the phase distribution as well as the phase shift amount of each interferogram. Okada [72] proposed an algorithm based on iteration to treat all the variables in the intensity distribution as ?unknowns?. The algorithm is very effective but it works only when the amount of phase shift error is small. When the shifting error becomes larger, a convergence problem occurs and accuracy decreases significantly. This is a rational for the proposed algorithm. 3.2.2.1 Enhanced iteration algorithm The algorithm is based on the non-linear least square method. An estimated intensity is defined for the cases of known phase values and known phase shifting amounts. An iteration process is performed until the sum of the squared difference between the estimated intensity and the measured intensity converges to a small value. 40 1. Step 1: Phase calculation under an assumption of no random noise The estimated intensity of the ith phase shifted image, eijI , can be defined as ( )eijijijjiijijjiijjiIABcosABcoscosBsinsin,i0,1,,,M1=+??+?=+?????????=? (3.33) where the subscript ?i? denotes the ith phase shifted image and the subscript ?j? denotes the individual pixel locations in each image. In the equation, ?i is the amount of phase shift of each frame, ?j is the unknown phase, Aij is the background intensity and Bij is the modulation amplitude. Assuming that there is no random noise in the patterns, Aij and Bij become single order tensors, i.e., A1j = A2j = AMj and B1j = B2j = BMj. Defining a new set of variables as jijaA= , jijjbBcos=?, jijjcBsin=??, the estimated intensity can be expressed as ijijj e ij sinccosbaI ??+??+= (3.34) In the classical least-square approach, an expression for the error Sj, accumulated from all the images, can be written as ( ) ( )?? ? = ? = ???+??+=?= 1M 0i 2 ijijijj 1M 0i 2 ij e ijj IsinccosbaIIS (3.35) In this step, it is assumed that i? have correct values. Then the least squares criteria require ?? ? ? ? ?? ? ? ? ? =?? =?? =?? 0cS 0bS 0aS j j j j j j (3.36) 41 The above equation yields ?? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ?= ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ???? ???? ?? ? ? ? ??? ??? ?? ? = ? = ? = ? = ? = ? = ? = ? = ? = ? = ? = 1M 0i iij 1M 0i iij 1M 0i ij j j j 1M 0i i 2 1M 0i ii 1M 0i i 1M 0i ii 1M 0i i 2 1M 0i i 1M 0i i 1M 0i i sinI cosI I c b a sincossinsin sincoscoscos sincosM (3.37) From above matrix equations, the constants aj, bj, cj can be determined and the unknown phase can be determined from ??? ? ??? ? ?=? ? j j1 j b ctan (3.38) 2. Step 2: Phase shift calculation under an assumption of uniform background In the second part of iteration, it is assumed that the background and modulation are constant for each frame. Then Aij and Bij become single order tensors, i.e., Ai1 = Ai2 = AiN and Bi1 = Bi2 = BiN. Defining another set of variables for each frame as iijaA? = , iijibBcos? =?, iijicBsin? =??, the estimated intensity can be expressed as jijii e ij sinccosbaI ???+???+?= (3.39) An expression for the error iS?, accumulated from all the pixels in the ith image, can be expressed as ( ) ( )?? ? = ? = ????+???+?=?=? 1N 0j 2 ijjijii 1N 0j 2 ij e iji IsinccosbaIIS (3.40) In this step, it is assumed that j? computed in Step 1 is correct. Then, the least squares criteria require 42 ? ? ? ? ?? ? ? ? =?? ?? =?? ?? =?? ?? 0cS 0bS 0aS i i i i i i (3.41) The result yields ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?= ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ???? ???? ?? ? ? ? ??? ??? ?? ? = ? = ? = ? = ? = ? = ? = ? = ? = ? = ? = 1N 0j jij 1N 0j jij 1N 0j ij i i i 1N 0j j 2 1N 0j jj 1N 0j j 1N 0j jj 1N 0j j 2 1N 0j j 1N 0j j 1N 0j j sinI cosI I c b a sincossinsin sincoscoscos sincosN (3.42) The constants ia? , ib? , ic? can be determined from the above equation. Then, the amount of phase shift in each frame can be determined from ??? ? ??? ? ? ??=? ? i i1 i b ctan (3.43) 3. Step 3 Repeat the steps until the phase shift values converge, i.e., a condition ?1, then 2k T 2k 1k T 1k k zr zr ?? ??=? 1kk1kk pzp ?? ?+= 48 (6) One scalar and two vector updates are performed: k T k 1k T 1k k Qpp zr ??=? kk1kk p?+?=? ? kk1kk Qprr ??= ? (7) If maxkk ? , or 0k rr ?< , the calculation stopped. Otherwise, go to Step (2). Surveys and comparisons indicate that PCG least-square iteration always gives the best performance for phase unwrapping. The disadvantage of PCG iteration method is that the unwrapping time is much longer than the above two methods. The details of the algorithm can be found in refs. [95] and [96]. 3.3.4 About the phase unwrapping methods The above three unwrapping algorithms are three typical methods among dozens of existing phase unwrapping methods. The sequential filling is a very fast method; it can be used for the fringe analysis of good-quality fringe patterns (e.g., interferometric moir? fringes). The minimum spanning tree method is much slower than sequential filling but generally yields reliable results for most of the photomechanics fringe analysis. The PCG least-squares iteration is a relatively slow method but it offers the best performance. Figure 3.9 is an example of phase unwrapping for a complicated fringe pattern using the three algorithms. It is evident that the PCG method yields the best result because the unwrapped phase is consistent with the original wrapped phase. 49 Some other phase unwrapping methods, such as Flynn?s minimum discontinuity method, L0 method and cellular-automata method, also offer good performance for phase unwrapping of fringe patterns. However, these methods are less practical and are not discussed here [92]. (a) Wrapped phase (b) Unwrapping using sequential filling (c) Unwrapping using minimum spanning tree (d) Unwrapping using PCG iteration Figure 3.9 Unwrapping of complicated phase using different algorithms 50 3.4 Calculation of fringe order gradient When an in-plane displacement field is studied, it is often required to obtain its derivative (or strain) to complete the deformation analysis. The strain distribution can be calculated from the full-field displacement field by U(x,y)N(x,y)k(x,y)(x,y)k xx2x ?????=?= ?pi? (3.51) where ? is an engineering strain component, U is the displacement field, k is a contour interval of the fringe pattern, ?N is the change of fringe orders in the fringe pattern and ?x is any gage length across which ?N is determined. In practice, the displacement field obtained from the experiment contains optical and electrical noise. The noise will not affect the displacement field much; however, it can result in large errors in gradient or strain calculation. For this reason, low-pass filtering or smoothing can be applied to the displacement fields to reduce the noise using general filtering techniques such as averaging filtering or median filtering. The displacement smoothing can also be implemented through a surface polynomial fitting. Similar to smoothing displacement fields before gradient calculation, an alternative technique is using the gage lengths of many pixels in the gradient or strain calculation. However, this arbitrariness in smoothing filer or gage length selection makes the accurate strain calculation difficult. Figure 3.10 illustrates an example of strain calculation using various gage lengths. The fringe pattern represents a horizontal direction displacement of an infinite plate with a hole subjected to a uniform tension. 51 As can be seen from the plots, a larger gage length smoothes out the strain plot but they mask the strain concentration at the hole boundary. A similar situation happens at the interface of dissimilar materials, where a strain gradient is extremely large. (a) Original images (b) Wrapped phase (c) Displacement 0 10 20 30 40 50 60 70 150 200 250 300 350 400 x position epsilon-xx (micro-strain) gage length: 1 gage length: 5 gage length: 20 Theoretical (d) Strain ?x plot along vertical center line Figure 3.10 Examples of strain calculation using different gage lengths 52 3.5 Limitations of automatic analysis The automatic fringe analysis techniques are very powerful tools to determine a full field phase (or fractional fringe order) map. The results can be used to plot a 2-D or 3-D deformed configuration. However, the automated analysis should be restricted to regions of singular material properties when a gradient is sought; they should not cross boundaries between materials of different properties because stresses (forces) are continuous functions across boundaries, so strains and displacement derivatives must be discontinuous wherever abrupt changes of properties (e.g., the elastic modulus) occur. Otherwise, the filtering or smoothing misrepresents the data near discontinuities. A similar argument has been made for a region with a stress concentration. This is the motivation of development of a semi-automatic process for strain calculation. 53 4. FRINGE ANALYSIS II: SEMI-AUTOMATIC ANALYSIS This chapter is devoted to semi-automatic processing techniques for fringe patterns. The exiting techniques are reviewed and the limitations are discussed. New hybrid algorithms are proposed and the preliminary results from the proposed algorithms are presented. 4.1 Single image: fringe centering method Figure 4.1 Flow chart of fringe centering method Image input Image pre-processing Start End Right ? Yes No Fringe centerline detection Fringe modification Boundary determination Fringe order determination Fringe order interpolation Experiment parameter calculation Results output 54 Fringe centering technique, also known as ?fringe skeletonizing?, was developed for the traditional manual fringe analysis, in which the fringe spacings are used to calculate displacements or strains. Before the phase measurement techniques became available, the fringe centering technique was the only processing tool available for the automatic analysis of interferograms. This class of techniques remains a vital element in the repertoire of fringe analysis methods. The fringe centering method is the only viable automatic fringe analysis technique for interferograms if only photographic records of the interferograms are available or the experimental event is dynamic such as impact testing. The flow chart of the fringe centering method is shown in Figure 4.1. The details of each step are explained below. 4.1.1 Fringe centerline detection In this method, only the integer fringe orders along the fringe centerlines are sought. The process of fringe centerline detection is to find points representing the fringe centerlines by eliminating all other parts of the fringes. There are two different methods for extracting the fringe centerlines from the fringe patterns [8]~[12]. The first one involves binarizing the fringe patterns and then skeletonizing the binary fringes; the other method detects the local maxima or minima of fringe intensities in a gray-scale. The peak detection methods are more sensitive to noise than their binary counterparts, but they offer the prospect of higher resolution detection of the fringe centers. 55 4.1.1.1 Binarization In the binarization process, the gray levels above or below a threshold value are truncated to the maximum or zero intensity, respectively, to convert the image into a binary intensity image [13]~[15]. The mean intensity of the fringe pattern can be chosen as the binarization threshold. When the image has an uneven background, the image can be segmented into small blocks and binarizing operation is applied to each block. After binarization, the geometric centerline of the black fringe is regarded as the fringe centerline. This assumption can give rise to a large error if the fringe intensity is asymmetric. Because of this disadvantage, the fringe binary method is not widely used. 4.1.1.2 Peak detection Peak detection [16]~[24] means finding the local maxima or minima of the gray-scale images. In this method, the whole image is subjected to a peak detection matrix and the gray level of the pixels in the matrix is reduced to zero value if a peak does not exist and to logical ?1? if a peak is found. Among many fringe centerline detection methods, a 5?5 window pixel peak detection scheme is one of the easiest and most effective methods. Figure 4.2 Pixel matrix of 5?5 and directions for fringe peak detection P-2 2 P-1 2 P0 2 P1 2 P2 2 P-2 1 P-1 1 P0 1 P1 1 P2 1 P-2 0 P-1 0 P0 0 P1 0 P2 0 P-2 -1 P-1 -1 P0 -1 P1 -1 P2 -1 P-2 -2 P-1 -2 P0 -2 P1 -2 P2 -2 X Y XY -XY 56 The fringe peak detection uses two-dimensional peak detection, locally performed within a 5?5 pixel matrix, as shown in Figure 4.2. With respect to the four directions shown in the figure, the peak conditions are defined as follows: For the X-direction, 1,21,20,21,01,00,0 PPPPPP ????? ++>++ (4.1) 1,21,20,21,01,00,0 PPPPPP ++>++ ?? (4.2) For the Y-direction, 2,12,12,00,10,10,0 PPPPPP ????? ++>++ (4.3) 2,12,12,00,10,10,0 PPPPPP ++>++ ?? (4.4) For the XY-direction, 2,11,22,21,11,10,0 PPPPPP ????? ++>++ (4.5) 1,22,12,21,11,10,0 PPPPPP ????? ++>++ (4.6) For the -XY-direction, 2,11,22,21,11,10,0 PPPPPP ???????? ++>++ (4.7) 1,22,12,21,11,10,0 PPPPPP ++>++ ?? (4.8) When the peak conditions are satisfied for any of two or more directions, the object point is recognized as a point on a fringe skeleton. 57 4.1.1.3 Fringe thinning In many cases, the fringe obtained from the image binary or the peak detection schemes is usually wider than the width of one pixel. Binary image thinning algorithms [28]~[30] are required to further reduce the fringe width and thus to obtain a true fringe centerline. Among many proposed algorithms, Rosenfeld thinning algorithm and Hilditch thinning algorithm are employed in this study. 4.1.1.4 Fringe centerline improvement Fringe patterns uniquely define the displacement fields. Because of the noise, the fringe centerlines obtained from above steps inevitably have undesired defects, such as broken fringe centerlines, cross-connected fringe centerlines and fringe centerlines with short branches. These defects must be eliminated before next processing. A series of corresponding algorithms to cope with these fringe centerline defects are utilized in this dissertation; these algorithms include automatic broken-line connection, automatic cross- connection line seperation, short fringe branch elimination, manual fringe connection and elimination, and so on [31]. Figure 4.3 shows an example of fringe centerline detection using binarization and peak detection methods. (a) is the original experiment horizontal field fringe pattern of a diametrical compression circular disc with a hole in the center; (b) is the fringe pattern after low-pass filtering of (a); (c) is the image after binarizing; (d) is the fringe centerline after thinning and (e) is the image after improving fringe centerlines; (f) is the overlap of original fringe pattern (a) and the fringe centerline pattern (e) where the color of the 58 centerlines were changed to be white. It is evident that the fringe centerlines are well recognized. Similarly, (g) through (j) are the results of using fringe peak detection method. (a) Original fringe pattern (b) Low-pass filtering of (a) (c) Binarization of (b) (d) Thinning of (c) 59 (e) Fringe improvement of (d) (f) Overlap of (a) and (e) (g) Peak detection of (b) (h) Thinning of (g) (i) Fringe improvement of (h) (j) Overlap of (a) and (i) Figure 4.3 Example of fringe centerline detection 60 4.1.2 Semi-automatic fringe orders assignment Fringe centerlines represent the contours of equal displacements; the fringe centerline itself does not contain information of fringe orders. Therefore, fringe order assignment is necessary for the whole-field fractional fringe order calculation [11][12]. The adjacent fringe orders differ by -1 or 1 except in zones of local maxima and minima. Based on this feature, a semi-automatic process to assign fringe orders has been developed to achieve fringe ordering. This semi-automatic process requires a simple human-computer interaction and the increasing or decreasing directions of some fringe orders should be known. The later requirement can be provided through a judgment based on the mechanical behavior of the testing specimen or the moving orientation of the fringes when small phase change is added into the experiment system. Figure 4.4 Example of semi-automatically fringe order assigning L1 L1? H1 H1? L2? L2 H3? H3 O H2 L3 L4 61 Figure 4.4 shows an example of the process. Suppose the order of a fringe located at L1 is 0 and the direction from L1 to H1 is known (increasing or decreasing). Then the orders of the fringes spanned from L1 to H1 can be determined. Now, fringe orders at the fringes that contain H1? and H2 are known; according to this information, fringe orders from H2 to L2 can also be determined automatically. The same procedure can be used for other segments except the circular fringe located in the center whose fringe order determination requires a human judgment. The fringe order at L1? is compared with the fringe order at L1 to double check whether the fringe orders are assigned correctly. Finally, a constant fringe order can be added to all the fringes to reflect the real fringe orders. 4.1.3 Fractional fringe order calculation: fringe order interpolation In the fringe centering method, only the integral fringe orders along the fringe centerlines are determined. Interpolation is required to obtain fractional fringe orders at every pixel. The following sections describe the interpolation methods used in this study. 4.1.3.1 1-D interpolation The most widely used 1-D interpolation algorithm is cubic spline interpolation, and it has been successfully applied to fringe order interpolations [9][12]. Cubic spline interpolation requires the boundary conditions to be known, which is the most critical limitation of this algorithm since the boundaries are often regions of interest. For this reason, a segment-by-segment curve fitting interpolation algorithm is proposed to avoid 62 the requirements of boundary conditions. The proposed algorithm is based on a continuous differential of the first derivative (2nd order differential of displacement). The proposed algorithm is illustrated in Figure 4.5. For an arbitrary point X, the fringe order at point X can be calculated by: ?? ? ?? ??+ ?? ? ?? ??+ ?? ??? ? ?+ ?? ? ?? ??+ ?? ? ?? ??+ ?? ??? ? ?= )XX)(XX( )XX)(XX)(X(f )XX)(XX( )XX)(XX)(X(f )XX)(XX( )XX)(XX)(X(f XX XX )XX)(XX( )XX)(XX)(X(f )XX)(XX( )XX)(XX)(X(f )XX)(XX( )XX)(XX)(X(f XX XX)X(f 3424 324 4323 423 4232 432 23 2 2313 213 3212 312 3121 321 23 3 (4.9) where X1, X2, X3, X4 are the fringe centerline points and their fringe orders f( Xi ) are known. It is obvious from equation (4.9) that the derivative of the fringe order is continuous. Actually, the 2nd order differential of the fringe order is also continuous. Figure 4.5 1-D interpolation XX1 X2 X3 X4 X f(X4) f(X3) f(X2) f(X1) f(X) 63 4.1.3.2 Limitation of 1-D interpolation The main disadvantage of 1-D interpolation is that interpolations along x-direction and y-direction do not yield the exactly same results. This implies that 1-D interpolation is not sufficient to describe the full-field experimental parameters. The error can be expected when the derivatives are calculated along one direction using data obtained from interpolation along another direction (e.g., calculation of strain ?x in the x-direction using the U-field displacement data that were obtained from the interpolation along y- direction). 4.1.3.3 New approach: Improved 1-D interpolation Instead of using the fixed Cartesian coordinates, the directions normal to the fringe orientations can be used to improve the efficiency of 1-D interpolation. For example, the interpolation along the virtual-lines in Figure 4.6 should yield better results for the points along those lines simply because fringe orders change more rapidly along them. Once the fractional fringe orders are obtained at the points of the virtual lines, the virtual-lines can be regarded as regular fringe centerlines. These lines are called as ?assistant fringes?; although they are not real fringe centerlines, they are useful for the fringe order interpolation. Figure 4.7 shows an example of 1-D interpolation. The fringe pattern represents a vertical displacement of a plate with a hole subjected to a uniform tension. The results clearly indicate effective performance of the improved 1-D interpolation. 64 Figure 4.6 Improved 1-D interpolation (a) (b) (c) (d) (e) (a) V-field fringe; (b)Fringe centerlines; (c) X-direction interpolation; (d) Y-direction interpolation; (e)Improved X-direction interpolation Figure 4.7 Example of improved 1-D interpolation 65 4.1.3.4 Proposed 2-D interpolation approach: tile-by-tile interpolation Since fringe pattern is a 2-D image, a 2-D fringe order interpolation seems more attractive. For simple and uniform fringe patterns, a global 2-D polynomial interpolation can be used. However, for the real engineering problems, the 2-D interpolation should be based on tile-by-tile processing. To use a proper number of fringe centerlines, the size of the tile should be adjusted according to the fringe density. After the 2-D polynomial interpolation within every tile is conducted, the displacement data at the edge of adjacent tile are compared and smoothly connected. Figure 4.8 shows the result of the fringes in Figure 4.7, obtained by using the 2-D interpolation. Currently, 2-D interpolation is still in somewhat immature state, more development efforts are need to complete the 2-D interpolation. Figure 4.8 Example of 2-D interpolation 66 4.1.4 Fringe order gradient or strain calculation Fringe order gradient or strain can be calculated using the interpolation function, namely differentiating the polynomial function used in interpolation. Figure 4.9 shows an U-field example of a tensile specimen with a hole. The theoretical stress (also strain) concentration factor at the top point of the hole is 3.0, and the fringe centering method gives 3.08. (a) Fringe pattern (b) Strain distribution Figure 4.9 Example of strain calculation 4.2 Overall limitations of fringe centering method Compared with automatic fringe analysis methods (such as Fourier transform and phase shifting techniques), the fringe centering method uses only the fringe centerlines 67 for the processing. Because interpolation is employed to determine the fractional fringe orders, the strain calculation is simple with the fringe centering method. Figure 4.10 illustrates an example of using the fringe centering method to process the thermally induced vertical displacement fringe pattern of a solder ball interconnection in an electronic packaging component; 1-D interpolation is used in this example. (a) Original Fringe pattern (b) Image filtering (c) Fringe centerlines detection (d) Fringe centerlines improvement (e) Fringe orders assignment (f) Displacement field 68 (g) Strain field Figure 4.10 Example of fringe analysis using fringe centering method Despite its advantages, applications of the fringe centering method are cumbered with some limitations. They include: (1) The fringe centerline location should be accurately determined; otherwise, the strain calculation might produce a significant error. Unfortunately, accurate locating of the fringe centerlines is often difficult for the fringe patterns with noise. (2) Semi-automatic fringe order assignment requires some knowledge on fringe orders. With a complex loading and/or geometry, fringe orders can be ambiguous. (3) The interpolation is implemented based on the fringe centerlines; however, in many cases, the number of fringe centerlines available for accurate interpolation is often insufficient. (4) For many engineering problems, the most important areas are near the geometrical boundaries. However, the fringe centerlines are not usually 69 available at the boundaries or very close to the boundaries. The gradient at the boundaries can be erroneous depending upon a choice of interpolation function. To improve the accuracy of fringe centerline locations, sub-pixel operations through magnifying the original fringe image can be adopted. This involves an interpolation of the image. Another effective solution to the above limitations is using optical/digital fringe multiplication (O/DFM) technique [25]. The O/DFM method uses a series of phase- shifted fringe images. The method provides multiplied fringe centerlines with high accuracy. An added benefit is its ability of automatic fringe ordering. 4.3 New hybrid O/DFM fringe centering method for multiple images A hybrid method combining the optical/digital fringe multiplication (O/DFM) method [25] and the fringe centerline method is proposed to cope with the limitations of the fringe centering method. 4.3.1 Optical/digital Fringe multiplication (O/DFM) In O/DFM, a series of n shifted patterns are utilized [25]. These shifted patterns are sharpened and combined into a single contour map, which exhibits n time as many fringes as the original pattern. Figure 4.11 shows a schematic illustration of O/DFM for n = 2. The O/DFM algorithm subtracts the intensities at each point, as illustrated in (b). Then it inverts the negative portions by talking absolute values, as illustrated in (c). The algorithm proceeds by truncating the data near 0Ir = and binarizing by assigning intensities of zero and one 70 to points below and above the truncation value, respectively. The result is graphed in (d). The result is a sharpened contour map that has twice as many contour lines as the number of fringes in the initial pattern. The sharpened contours occur at the crossing points of the complementary graphs, where the intensities of the complementary patterns are equal. Note that any noise or other factor that affects the two patterns equally has no influence on the locations of the crossing points. (e) illustrates a fringe multiplication of 6. Figure 4.12 shows an example of the O/DFM method for n = 4 on a whole field basis, where the fringes represent the same fringe pattern in Figure 4.10. Compared with the conventional fringe centerline detection shown in Figure 4.10 (c), the O/DFM offers much higher quality of the fringe centerlines (Figure 4.12 d) and fringe centerline improvement is not required. In addition, the fringe orders can be assigned automatically. Figure 4.11 Schematic illustration of O/DFM 71 (a) The initial fringe patterns (b) O/DFM subtraction pattern (c) Sharpened contours (d) Fringe centerlines (e) Displacement field (f) Strain field Figure 4.12 Example of O/DFM 72 4.3.2 Hybrid approach The O/DFM method overcomes several limitations of the fringe centering method. The combination of the O/DFM and the fringe centering method makes the semi- automatic fringe analysis more effective. In the hybrid approach, the O/DFM method is first employed to obtain the accurate high-quality fringe centerlines. Then the fringe interpolation proceeds to obtain fractional fringe orders at every pixel. Besides noise reduction, the fringe orders can be determined automatically using the phase shifting algorithm described before. The displacement interpolation can be performed with higher fidelity because of the enhanced displacement sensitivity (i.e., more fringe contours). Figure 4.13 shows the results of the fringes in Figure 4.9, processed by the hybrid O/DFM fringe centering method. The strain concentration factor of 2.98 was achieved. (a) Fringe patterns (b) Strain distribution Figure 4.13 Example of strain calculation 73 4.3.3 Advantages and disadvantages of hybrid semi-automatic analysis 4.3.3.1 Comparisons with automatic analysis method Compared with automatic fringe analysis methods (i.e., Fourier transform method and phase shifting method), the hybrid semi-automatic method has the following advantages: (1) The hybrid approach does not require that the intensities of fringe patterns should be cosinusoidal. (2) The hybrid approach uses fringe order interpolation for the displacement and strain calculations; the effect of the noise is much smaller. (3) The hybrid approach is suitable for problems containing geometrical as well as material discontinuities. The only disadvantage of hybrid approach is that interpolation algorithm is required. 4.3.3.2 Comparisons with conventional fringe centering analysis method Compared with the conventional fringe centering method, the hybrid semi-automatic approach has the following advantages: (1) Fringe centerline detection is less sensitive to the noise; the fringe centerline location is determined more accurately. (2) The number of fringes is multiplied; this increases accuracy in fringe order interpolation. (3) Fringe orders can be assigned automatically. The disadvantage of the hybrid method is that it requires more than one image. 74 5. SOFTWARE DEVELOPMENT Software design and development are important portions of this dissertation. The software is developed using Microsoft Visual C++ 6.0 and Microsoft Visual Studio .Net. The software is based on Windows Graphic User Interface (GUI) and can be executed on Microsoft Windows 9X/NT/2K/XP operating systems. Figure 5.1 shows the main architecture of the fringe analysis software schematically. The software includes the general digital image processing algorithms and all the algorithms discussed in this dissertation. The main features of the software include but not limited to: (1) General features ? Create, open, close, save image or images ? Image preview ? Print and print preview ? Undo and redo ? Cut, copy and paste ? Status and tool bars ? Display setup (2) General Image processing ? Adjust color and balance color ? Process multiple images simultaneously 75 ? General drawing tools: line, ellipse, rectangular, eraser, etc. ? Resize image ? Crop image ? Rotate image ? Zoom in and out ? Create user-defined color palette ? Create and copy boundary (3) Pre-processing: image filtering ? Spatial neighborhood-averaging filter ? Spatial median filter ? FFT and DFT low-pass filters ? FFT and DFT high-pass filters ? Self-adaptive filters (4) Automatic fringe analysis ? General Fourier transform ? Carrier Fourier transform ? Conventional phase-shifting algorithms (automatic frame detection) ? Random phase-shifting algorithm ? Four kinds of advanced phase-unwrapping algorithms (sequential filling, minimum spanning tree, PCG least-squares iteration and cellular-automata) ? Displacement field smoothing ? Strain calculation with different differential gaps ? Shear strain calculation 76 (5) Semi-automatic fringe analysis ? Image binarization ? Fringe peak detection ? O/DFM and O/DFM2 ? Automatic fringe thinning ? Automatic and semi-automatic fringe improvement techniques ? Automatic fringe searching ? Semi-automatic and automatic fringe order assignment ? 1-D cubic spline interpolation ? 1-D segment-by-segment curve fitting interpolation ? Improved 1-D interpolation: assistant fringes selection ? Global 2-D interpolation ? Tile-by-tile 2-D interpolation ? Strain calculation with refinement ? Shear strain calculation (6) Post-processing: display and output ? Display in gray-scale or color scale ? Section display ? Isoline display ? 2-D in-plane deformation display ? 3-D display (rotation and hiding) ? Data and image compression and decompression ? Data and image storage 77 ? Data interface with other software ? Image and result printing ? Report creation (7) Other features ? On-line Help ? Tip of the day ? Operation hint ? Fringes simulations (more than 10 kinds of mechanics experiments) ? Pixel value plotting ? Software demo ? Standard windows application software installation The software is menu-driven and mouse-driven. Figure 5.2 illustrates some screen captures of the software. 78 Figure 5.1 Schematic architecture of the fringe analysis methods in software Input Single or Multiple image(s)? Fringe filtering Fringe centering Fourier transform ODFM hybrid Phase shifting Fringe center detection Multiple-fringe center detection Phase calculation Phase calculation Specification of the boundary Fringe centerline thinning Fringe centerline improvement Fringe order assignment Displaceme Strain Phase unwrapping Fringe order interpolation Save and display S M 79 (a) Strain field calculation (b) Multiple images processing and operation hint 80 (c) 2-D deformation calculation (d) 3-D warpage calculation Figure 5.2 Examples of screen captures of the software 81 6. APPLICATION I: OUT-OF-PLANE SHAPE AND WARPAGE MEASUREMENT Applications to shape determination and warpage measurement are presented. Different experiment methods and fringe processing techniques are employed based on the specific requirements for each measurement. Among these applications, a new infrared diffraction interferometer is proposed for out-of-plane co-planarity measurement of high density solder bump pattern. 6.1 Infrared diffraction interferometer for co-planarity measurement of high density solder bump pattern An IR diffraction interferometer is proposed for co-planarity measurement of high- density solder bump patterns. The method utilizes long wavelength (? = 10.6 ?m), coherent infrared laser light, which serves to reduce the apparent roughness of test objects, and enables the regularly spaced solder bump arrays to produce well-defined diffracted wavefronts. A single diffracted wavefront is isolated by an optical system and directed to interfere with a reference wavefront to produce a whole-field map of bump topography. An optical configuration similar to classical Fizeau interferometry is implemented to prove the concept. Digital moir? fringe processing techniques are used in the experiment data analysis. 82 6.1.1 Introduction The ability to accurately monitor bump co-planarity is critical to the prevention of non-wetting of chip level, high density solder interconnects during assembly. Bump co- planarity has traditionally been measured through existing methods developed for topographical mapping. These include projection moir?, shadow moir? and white light surface profilometers [102]~[106]. The first two methods can map a relatively large area on a whole-field basis, but their application to high-density bump co-planarity has been limited due to inherent limitations in spatial resolution. The white light surface profilometer is a scanning, point-measurement technique that provides high-resolution measurements of surface topography. With sophisticated algorithms and automation techniques, the point-wise profilometer data can be stitched together to generate a map of surface topography. However, this method does not naturally generate a quick turn and real time map of surface topography as provided by whole-field measurement techniques. This research work proposes a new method to cope with the above limitations; namely, a whole-field interferometric method to measure the co-planarity of high-density bump patterns. The method utilizes coherent, infrared (IR) laser light source of long wavelength (?= 10.6 ?m) (1) to relax the surface roughness limitations of visible light interferometers, and (2) to enable the regularly spaced solder bump arrays to produce well-defined diffracted wavefronts. In the proposed scheme, a single diffracted wavefront is isolated and it is directed to interfere with a reference wavefront to produce a whole-field map of bump topography. The principles of the method and selected examples are described. 83 6.1.1.1 Shifting to a longer wavelength By shifting to a longer wavelength, a surface, which typically acts to diffusely scatter a collimated visible light beam, may become specular. This has been the precise motivation for developing long IR light interferometric systems to measure the topography of non-specular surfaces [107]~[112]. Based on a bi-directional reflectance distribution function (BRDF), the specular component dIspec, which is responsible for the sharp mirror-like reflection generated by a Gaussian, isotropically rough surface, can be expressed as [113]~[115] 2 g specidIFeSI ?= (6.1) where F represents the Fresnel reflectivity, S is a geometrical shadowing function and Ii is the incident radiance. The function g, which depends on the rms surface roughness ?, is given by 24cos g pi????= ????? (6.2) where ? is the angle of incidence, and ? is the wavelength of light. According to equation (6.1), the specular component increases as the value of g decreases. The influence of surface roughness on specularity is then captured by the ?apparent? roughness factor: r 4cospi???= ? (6.3) The factor, ? cos?, represents the projection of the rms surface roughness into the incident light direction. 84 A surface therefore appears rough or smooth to an incident light wave depending upon whether the apparent roughness is large or small. It is interesting to note that all surfaces appear smooth as the incident light approaches grazing angles, that is, as ? approaches 90?. Alternatively, for a given degree of surface roughness and a fixed angle of incidence, specular reflection can be increased by lowering the ratio ?/?, i.e., by an increase in the wavelength of the light source. In the proposed approach, a ? = 10.6 ?m light beam, generated by a CO2 laser, reduces the apparent surface roughness of the test sample by a factor of 20 at normal incidence (? = 0?), compared with a wavelength in the middle of the visible spectrum (green light with 0.5 ?m). The surfaces defined by the solder bump array and underlying substrate, appear optically rough to visible light, yet produce specular reflections when illuminated by the longer wavelength IR radiation. In addition, the long wavelength light can interact with the regularly spaced arrays to produce well-defined diffracted wavefronts. These two features, namely the conversion of a diffusive surface to a specular surface and the diffraction phenomenon associated with it, form the basis of the proposed method. 6.1.1.2 Problems associated with classical IR interferometry Traditional IR interferometry in the literature has been implemented using the optical configurations of classical two-beam interferometry; namely, Twyman/Green interferometry [108]~[110] and Fizeau interferometry [111][112]. In both configurations (Figure 6.1), a beam splitter directs one-half of the incident light onto the specimen surface along a direction virtually normal to its surface, while the other half is reflected 85 by an optically flat reference surface. The wavefront from the specimen (?2) interferes with the wavefront from the reference mirror (?1). The resulting interference pattern is then viewed upon the surface of interest to produce a contour map of the z coordinate of the specimen surface. In Fizeau interferometry (Figure 6.1 b), the reference path and the active path are identical, whereas the active path is perpendicular to the reference path in T/G interferometry (Figure 6.1 a). Consequently, the Fizeau interferometer is much easier to tune. This feature is especially advantageous to an IR optical system inasmuch as the IR light is not visible during normal operation [111][112]. With the configuration shown in Figure 6.1 b, fringes are usually visible on the monitor with small or no adjustment after the specimen is aligned to be parallel to the optical flat through a simple visual inspection. (a) Twyman/Green interferometry 86 (b) Fizeau interferometry Figure 6.1 Optical configuration of classical interferometry When the substrate and high-density solder bump pattern are normally illuminated by the long IR light, a complex light field is generated. The situation can be best described by the combined effect of two separate sets of uniform structures, as illustrated schematically in Figure 6.2; the solder bump array and the uniform spacing of the exposed underlying substrate. These two periodic structures behave like separate diffraction gratings since each surface has different reflectivity; a phase-type grating (solder bump; referred to as bump grating) and an amplitude-type grating (uniform 87 spacing; referred to as substrate grating). The wavefront profile of the specular 0th order diffraction is thus affected by the reflectance and topography of both arrays. Figure 6.2 Schematic illustration of reflection from the bump arrays and the underlying spaces In general, the surface of the substrate is not flat and it is different from the top surface profile of the bump array. Thus, the 0th diffraction order from the substrate grating produces a warped wavefront, which interferes with the 0th diffraction order generated by the bump grating and produces a wavefront that is confounded and un- interpretable. Consequently, the optical configuration of classical two-beam interferometry requiring normal illumination and specular reflection cannot be used to map the top surface of the bump profile. The wavefront from the bump array must be isolated, and the following optical configuration is proposed to achieve the condition. 88 6.1.2 Infrared diffraction Fizeau interferometry 6.1.2.1 Basic concept When an amplitude grating is illuminated by a collimated beam, most of the diffracted radiation is concentrated along the direction of specular reflection, or the 0th order direction, from the surface normal [116]. Consequently, the influence of the substrate grating, which acts as an inefficient amplitude grating, will diminish drastically if a diffraction order other than the 0th order is viewed. Figure 6.3 Basic concept of IR diffraction interferometry The basic principle of the proposed method is illustrated in Figure 6.3. A collimated IR beam illuminates a partially reflective optical flat placed near the specimen. One side of the optical flat is coated with an anti-reflection coating to avoid ghost patterns. A 89 portion of the collimated beam is reflected from the uncoated surface of the optical flat and the transmitted beam is diffracted from the specimen surface. In the set up, the specimen is tilted in such a way that a particular diffraction order becomes parallel to the reference beam. The diffracted wavefront interferes with the reference beam to produce an interference pattern. From the grating equation with and with respect to proper sign conventions [117], msinsinmp ??=??+ (6.4) where ? is the tilt angle, m is the diffraction order, p is the pitch of the bump array and ?m is the angle of the mth diffraction order. The incident angle and the diffraction angle are measured with respect to the axis perpendicular to the plane of bump array. The condition that a diffracted wavefront becomes parallel to the reference wavefront can be achieved when the mth diffraction order has an angle of diffraction identical to the tilt angle; i.e., ?m = ?. From equation (6.4), the condition provides the tilt angle requirement as 1 msin 2p ? ????= ?? ?? (6.5) When the tilt angle satisfies the above condition, the mth order beam will now be directed back along the incident axis, which is coincident with the planar reference beam generated by the optical reference flat. The 0th order beam (containing the un-wanted substrate reflection) now makes an angle of 2? with respect to the incident axis and is therefore completely diverted from the interference path. 90 As the diffraction order increases, the intensity of the diffracted wavefront decreases rapidly. Consequently, the high orders are not desired in practice. Considering the first order diffraction, the required tilt angle becomes 1sin 2p ? ????= ?? ?? (6.5)? The required tilt angle for various bump pitches is plotted in Figure 6.4. The plot clearly shows that the shallow first order diffraction angles result across the domain of typical bump pitch values. Only a few degrees of specimen tilt angle are required in order to obtain a spatially filtered interference pattern. 0 5 10 15 20 25 30 35 0 50 100 150 200 250 300 350 400 Bump pitch (mm) Tilt angle ( ? ) Figure 6.4 Required tilt angle for different bump pitches when the first diffraction order is used 91 6.1.2.2 Governing equation The analysis used to determine the difference in optical path lengths follows the approach of Refs. [118]. Figure 6.5 illustrates details. Consider any point P on an idealized bump surface (i.e., perfectly uniform array with an equal bump height). If the bump surface deviates from the ideal position, point P moves to a new location P?, where W and U are the out-of-plane (planarity) and in-plane (irregular pitch) components of deviation, respectively. Figure 6.5 Changes of optical path length when point P moves to P? In the figure, a plane wavefront transmitted through the optical flat is incident upon the bump surface at ?. The optical path from the light source is the same for every ray in 92 the incident beam. The change of path length (?OPL) with respect to the incident beam, caused by the movement of the point, is analyzed in the figure, with the result OPL(x,y)2W(x,y)cos2U(x,y)sin?=?+? (6.6) where ?OPL represents the fringe order N at each point in the interference pattern by OPL(x,y)N(x,y) ?= ? . In a well manufactured bump array, variations in a bump pitch are usually negligible ( U0? ). Assuming for the moment that we can neglect the in-plane term, 2U(x,y)sin ? , the out-of-plane position, W, can be expressed as ( ) ( )Wx,yNx,y2cos?= ? (6.7) Equation (6.7) is identical to the governing equation of classical Fizeau interferometry with a small inclined angle of illumination ?. Using equation (6.5)?, the governing equation of the IR diffraction Fizeau interferometry can be written as ( ) ( )2Wx,yNx,y 21 2p ?= ???? ???? (6.8) The above equation defines a contour interval as 2 21 2p??????? ?? displacement per fringe order. For typical bump pitch values ranging from 150 to 500 ?m, 1p21 2 ??? ? ? ??? ? ?? and the contour interval remains virtually constant as 5.3m 2 ? =?. 93 The bump pitch was assumed uniform in the derivation of equation (6.8) was derived. In reality, the bump pitch may be slightly irregular due to inherent manufacturing variability. If the bump pitch is irregular, the shape of the diffracted wavefront (and thus the fringe order) is influenced by perturbations in the bump pitch as well as the bump surface topography. Recalling equation (6.6), the apparent fringe order N? can be defined as OPL(x,y)2W(x,y)cos2U(x,y)sinN(x,y) ??+?? == ?? (6.9) Using equation (6.5), the W position can be expressed as ( ) ( ) ( ) ( )22mWx,yNx,yUx,yNx,yp mm2121 2p2p ???=?=?? ?????????? ???????? (6.10) Therefore mN(x,y)N(x,y)U(x,y) p? ?= (6.11) Equation (6.11) defines an error in reading a fringe order, caused by the irregular bump pitch. The error increases linearly as the local bump pitch irregularity as well as the diffraction order increases. In a typical manufacturing environment, the bump pitch deviation of 0.5% is acceptable (e.g., 1.0 ?m deviation for 200 ?m pitch), which will disturb the fringe order only by 1/200, or a co-planarity error of less than 0.03 ?m when the first order diffraction (m = 1) is used. For this reason, equations (6.7) and (6.8) are regarded as the governing equations of the proposed technique. 94 6.1.2.3 Optical configuration The complete optical setup is illustrated in Figure 6.6. An air-cooled all-aluminum CO2 laser (Synrad: Model 48-I) was used as a coherent light source. It produced a highly stable TEM00 single mode. The laser emitted a vertically polarized light with the central output wavelength of 10.6 ?m. Figure 6.6 Optical configuration of IR diffraction interferometry The laser beam was first expanded by a plano-convex lens. A beam splitter directed the light to a collimating lens, which also served as a field lens. An optical flat was placed next to the collimating lens. The two interfering beams were collected by an imaging system. The imaging system consisted of an imaging lens and an IR CCD camera (Electrophysics: Model PV-320). All the optical elements described above were fabricated from ZnSe. 95 All surfaces except one side of the optical flat were coated with anti-reflection coating. ZnSe has a refractive index of 2.4208 at ? = 10.6 ?m, which yields reflectivity of about 17% at normal incidence. This high reflectivity of the uncoated surface is an added advantage, which can accommodate surfaces with a wide range of reflectivity while maintaining good fringe contrast [111]. An aperture plate was placed at the focal plane of the field lens (or collimating lens) to isolate the two interfering beams. It was accomplished by the following simple procedure: (Step 1) Adjust the optical flat and the specimen to make them parallel to each other (Step 2) Arrange the aperture so that the reference beam and the 0th order diffraction from the specimen pass through the aperture, and (Step 3) Tilt the specimen slowly with respect to the axis normal to its surface until a new fringe pattern is observed. Representative fringe patterns obtained from the above procedure are shown in Figure 6.7. The images in (a) and (b) represent an interference pattern obtained from the 0th order diffraction (Steps 1 and 2) and the first order diffraction (Step 3), respectively. The images clearly demonstrate the effect mentioned in the previous section; i.e., that in image (a) the wavefront from the bump array is combined with the wavefront from the substrate to produce a very complex and uninterpretable 0th order wavefront. 96 (a) (b) Figure 6.7 Fringe patterns obtained from (a) the conventional IR Fizeau interferometry and (b) the proposed method 6.1.3 Co-planarity of flip-chip package bump pattern It is a standard practice to manufacture high density solder bumps directly onto an organic substrate. The bumping process on the substrate is achieved through a low cost reflow process. The proposed method was implemented to measure co-planarity of a bump array on a flip-chip substrate. The cross-sectional view of the specimen is similar to the one shown in Figure 6.2 and the pitch of the array was 200 ?m. The fringe patterns produced by the proposed method are shown in Figure 6.8. A series of four fringe-shifted interferograms was recorded with sequential changes of phase by a constant increment of pi/2. This was accomplished by a successive movement of the optical flat in the direction normal to the specimen by 8cos? ? . A high precision motorized stage was employed to achieve the desired fringe shifting accuracy. 97 (a) 0 (b) pi/2 (c) pi (d) 3pi/2 Figure 6.8 Four phase-shifted images obtained from a flip-chip bump arrays on an organic substrate The four phase-shifted fringe patterns were analyzed by using digital fringe image random phase shifting technique introduced in previous chapter. The real phase shift amounts are detected to be 0, 93?, 179? and 283?, respectively; sequential filling 98 unwrapping algorithm is employed to get the full-filed phase map. The maps of the wrapped phase and unwrapped phase are shown in Figure 6.9 and the resultant 3-D map is shown in Figure 6.10. The deviations along the two lines are also plotted in Figure 6.11. The maximum deviation was seen in the line AA? and its magnitude was approximately 7 ?m. Considering a typical bump height of 60 ?m, it was a significant deviation, which must be known before assembly process. (a) wrapped phase (b) unwrapped phase Figure 6.9 Results of fringe processing 99 Figure 6.10 Three dimensional representation of bump co-planarity (a) 100 0 1 2 3 4 5 6 7 8 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Distance from the center along AA' (normalized by width) Displacement ( mm) (b) 0 1 2 3 4 5 6 7 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 Distance from the center along BB' (normalized by diagonal) Displacement ( mm) (c) Figure 6.11 Deviations along (b) AA? and (c) BB? 101 6.1.4 Summary An IR diffraction interferometer technique has been proposed for co-planarity measurement of high-density solder bump patterns. The method utilized coherent IR laser light of long wavelength (? = 10.6 ?m). The long wavelength served to reduce the apparent roughness of the solder bump and substrate surfaces, and thus increased specular reflection significantly. This roughness reduction enabled the high-density bump arrays to produce well-defined diffracted wavefronts. A single diffracted wavefront was isolated by an optical system and it was directed to interfere with a reference wavefront to produce a whole-field map of bump topography. An optical configuration similar to a classical Fizeau interferometer was implemented to successfully document the topographical map of a solder bump array. The method was proven effective for co-planarity measurement of high-density bump arrays. Digital fringe image processing techniques, specifically, phase shifting technique and random phase shifting technique, provide the robust tool to automatically obtain the out-of-plane shape ad warpage for the co-planarity measurement. 6.2 Warpage measurement of electronic packaging components Electronic packaging is progressing toward integrating more devices with more functions into a smaller confined space, while requiring a high level of reliability. New electronic components, materials, fabrication processes and configurations are emerging to achieve these goals. From the structural point of view, efforts are being made to stack as many layers of materials as possible within a confined space. These layers consist of 102 different materials, such as, Cu, Al, silicon, polymers, ceramics, and solder. The materials possess different mechanical, thermal, and hygroscopic properties. After combinations of these materials are laminated, each material will exhibit a unique behavior under the variation of mechanical, thermal, and hygroscopic loading. Warpage is a global effect of interfacial stress and displacement. Accumulated interfacial stress and relative displacement will eventually delaminate and fail the structure. Warpage can also mis- registration and non-contact between a component and its substrate. For high-density interconnection, such as flip-chip (FC), chip scale packaging (CSP), and ball grid array (BGA), maintaining a flat surface to achieve high solder connection yield is vitally important [120]. Photomechanics is taking a leadership role for the warpage measurement of electronic packaging components. The general experimental techniques for out-of-plane displacement measurement include Twyman-Green interferometry [3][121][122], Fizeau interferometry, far infrared Fizeau interferometry [111][112], shadow moir? [123]~[127] and projection moir? [120]. Among these methods, Twyman-Green interferometry and shadow moir? are the most widely used techniques; their typical sensitivities are 0.316 ?m and 25 ?m respectively. 6.2.1 Warpage measurement using Twyman-Green interferometry The optical configuration of Twyman-Green interferometry was introduced in Figure 6.1 (a). The interference pattern is a contour map of the z coordinate of the specimen surface, where the contour interval is half of the wavelength of the laser light. The governing equation is 103 ( ) ( )y,xN2y,xW ?= (6.12) where W is the out-of-plane displacement, N is the fringe order and ? is the wavelength of the laser light. Figure 6.12 shows an example of warpage measurement of a tape-automated bonding (TAB) package using Twyman-Green interferometry. In this experiment, a series of 10 frame phase shifted fringe images were captured randomly (a). The random phase shifting algorithm was employed to obtain the wrapped phase (b). After unwrapping using the sequential filling method, the unwrapped phase map was obtained (c). It is worth noting that the defects in the original fringe patterns (a) were masked, which can be seen from the wrapped phase map (b); these defects were automatically corrected using non-linear interpolations after the unwrapping. The 3-D warpage was shown in (d) (a) Fringe image patterns 104 (b) wrapped phase (c) unwrapped phase (d) 3-D warpage Figure 6.12 Warpage measurement of a TAB package 105 6.2.2 Warpage measurement using shadow moir? method Figure 6.13 illustrates the basic concepts of shadow moir? [117]. The specimen is prepared by spraying it with a matte white paint. A linear reference grating of pitch g is fixed adjacent to the surface. A light source illuminates the grating and specimen, and the observer or camera receives the light that is scattered in its direction by the matte specimen. The explanation of shadow moir? assumes that the shadow of the reference grating is also a grating; the superposition of the shadow grating and reference grating forms a moir? pattern. From the cross-sectional view of Figure 6.13, the governing equation of shadow moir? can be obtained as ( ) ( )y,xNtantan gy,xW ?+?= (6.13) where W is the out-of-plane displacement, N is the fringe order, g is the pitch of the reference grating, ? and ? are the angles of incoming and outgoing rays, respectively. The frequency of the grating used in shadow moir? is usually less than 40 lines/mm (corresponding pitch is 0.025 mm/line), thus its sensitivity is much lower than Twyman- Green interferometry. This relative low sensitivity makes shadow moir? for measurement of relatively large warpage. Figure 6.14 shows another example of 3-D surface measurement of a quarter coin. This measurement cannot be implemented by using Twyman-Green interferometry because its sensitivity is too high for this surface measurement. The random phase shifting and the PCG iterative phase unwrapping algorithm were employed for the fringe analyses in this example. 106 Figure 6.13 Principle of shadow moir? method (a) Fringe patterns (b) Wrapped phase 107 (c) Unwrapped phase (d) 3-D shape Figure 6.14 Surface measurement of a quarter coin Figure 6.15 shows an example of warpage measurement of a plastic ball grid array (PBGA) package using shadow moir?. In this experiment, a series of four frame phase shifted fringe images were captured (a). The random phase shifting algorithm was employed to obtain the wrapped phase (b). After unwrapping using the minimum spanning tree method, the unwrapped phase map was obtained (c). The 3-D warpage was shown in (d). (a) Fringe patterns (b) wrapped phase 108 (c) unwrapped phase map (d) 3-D warpage Figure 6.15 Warpage measurement of a PBGA package 6.2.3 Warpage measurement using far infrared Fizeau interferometry (FIFI) The optical configuration of Fizeau interferometry was introduced in Figure 6.1 (b). Similar to Twyman-Green interferometry, the reflection beam from the specimen and the reflection beam from the optical flat meet again and they are collected by the camera. The interference pattern seen by the camera is the contour map of separation between the warped wavefront ?2 and the plane wavefront ?1. The interference pattern is a contour map of the z coordinate of the specimen surface, the governing equation is ( ) ( )y,xNcos2y,xW ??= (6.14) where W is the out-of-plane warpage or displacement, N is the fringe order and ? is the wavelength of the laser light. By employing a far infrared light with a very long wavelength (10.6 ?m), the surface that are regarded as rough under visible light can be tested. 109 Figure 6.16 shows an example of warpage measurement of a flip-chip plastic ball grid array (FC-PBGA) package using far infrared Fizeau interferometry (FIFI) [6]. In this package, a square chip (12 mm? 12 mm) was mounted on a substrate (12 mm? 12 mm). The package was virtually flat at the underfill curing temperatures; however, the large mismatch in CTE caused the package to bend as the temperature decreased. The fringe patterns obtained at 150?C, 100?C and room temperature are shown in (a)-(c). The laser wavelength is 10.6 ?m, therefore the contour interval is 5.3 ?m per fringe order. The fringes were analyzed using the fringe centering method. The fringe centerline is shown in (d); the 2D and 3D warpage maps at room temperature are shown in (e) and (f), respectively. Figure 6.16 (f) reveals a significant bending in both chip and substrate. Irregularities shown in the substrate are not caused by optical noise. Instead, they represent the heterogeneous deformation of the multi-ply glass/epoxy composite. (a) 150?C (b) 100?C 110 (c) Room temperature (d) Fringe centerline (e) 2-D warpage (f) 3-D warpage Figure 6.16 Warpage measurement of a FC-PBGA package 111 7. APPLICATION II: IN-PLANE DISPLACEMENT AND STRAIN MEASUREMENT In this chapter, Applications of in-plane displacement and strain measurement are presented. Among these applications, an inverse method to determine elastic constants is a new method and it is described in details. 7.1 Inverse method to determine elastic constants using circular disc and moir? interferometry An inverse method using a circular disc in diametrical compression is proposed and used for simultaneous determination of two elastic constants, E and n, from a single displacement map. Moir? interferometry combined with the phase-shifting technique provides a full-field displacement field. An over-deterministic approach using the nonlinear least squares method is implemented to fit the experimentally determined displacements to the theoretical solution. An implementation guideline is provided considering the effects of accidental rigid-body motions, random noise and imperfect position of the origin. Accuracy of the proposed method is verified experimentally. 7.1.1 Introduction A circular disc in diametrical compression is an experimental configuration that is easy to machine and load. With the well-established theoretical displacement fields 112 [128]~[132], a classical coefficient inverse approach can be implemented to determine material constants and/or applied load from the experimentally determined displacement fields [133][134]. Moir? interferometry measures in-plane displacements, U and V. The method is characterized by a list of excellent qualities, including full-field technique, high measurement sensitivity, high spatial resolution and high signal-to-noise ratio [135]. The data are received as whole-field interference fringe patterns, or contour maps, of the displacement fields. Because of the high sensitivity and abundance of data, reliable strain distributions?normal strains and shear strains?can be extracted from the patterns. For most practical applications, where the analyses are designed to investigate specific characteristics of the structure, quantitative results are extracted only at designated locations ? at points, or along lines in the specimen ? and whole-field analysis techniques are not employed. In the original attempts of the inverse approach using moir? fringes [133][134], only a limited number of displacement data were utilized for the analyses. The full?field displacement information had not been utilized completely. The full-field data are used advantageously in the current work In this dissertation, an inverse approach is proposed to determine two elastic constants from the whole-field displacement information provided by a moir? pattern. With the aid of a digital image processing scheme, displacement information is obtained at every point, which allows use of an over-deterministic analysis by the nonlinear least square method. Young?s modulus and Poisson?s ratio are determined simultaneously from a single moir? fringe pattern. 113 7.1.2 Background 7.1.2.1 Moir? interferometry with phase shifting Moir? interferometry has matured rapidly to emerge as an invaluable tool, utilized by many industrial and scientific applications [135][136]. In the method, a high-frequency cross-line grating on the specimen, initially of frequency fs, deforms together with the specimen. Two mutually coherent beams of laser light form a virtual reference grating, which interacts with the deformed specimen grating to produce a moir? fringe pattern. The principle of moir? interferometry is depicted schematically in Figure 7.1. The moir? patterns are contour maps of the U or V displacement fields, i.e., the displacements in the x and y directions, respectively, of each point in the specimen grating. The relationships, for every x,y point in the field of view, are x s y s 1U(x,y)N(x,y) 2f 1V(x,y)N(x,y) 2f = = (7.1) In routine practice of moir? interferometry, fs = 1200 lines/mm. In the fringe patterns, the contour interval is 1/2fs, which is 0.417 ?m displacement per fringe order. 114 Figure 7.1 Schematic diagram of moir? interferometry As stated in previous chapter, the phase shifting method utilizes a series of phase- shifted digital images to determine the phase at every point in the moir? pattern. Recalling it here, with succesive phase shifts of ? = pi/2, one can obtain the intensity distributions of four phase-shifted images as ( ) ( ) ( ) ( )imaIx,yIx,yIx,ycosx,yi,i0,1,2,32??pi??=+?+=?????? ?? (7.2) where Ii is the intensity at a point (or pixel) in the moir? pattern, Im is the mean intensity, Ia is the intensity modulation amplitude, ? is the angular phase information of the moir? 115 pattern, and (x,y) represents all the points in the x-y plane of the object and the moir? pattern. Then, using the well-established phase shifting algorithm introduced previously, the phase can be expressed in a simple form as ( ) ( ) ( )( ) ( )? ? ? ?? ? ? ?= y,xIy,xI y,xIy,xIarctany,x 20 13f (7.3) After the unwrapping process, the phase ? represents the fringe order N at each point of the pattern by ( ) ( )x,yNx,y 2?= pi . Then the whole-field displacement can be obtained with equation (7.1). 7.1.2.2 Displacement fields of circular disc in diametrical compression A circular disc in diametrical compression is shown schematically in Figure 7.2. Using the approach of Timoshenko and Goodier [138], the in-plane displacement components for the plane stress condition can be expressed as xxx yyy 1uAB EE 1uAB EE ?=+ ?=+ (7.4) where ux and uy : theoretical in-plane displacement components, E : Young?s modulus, ? : Poisson?s ratio, ( ) ( ) ( ) ? ? ? ?? ? ??????+? pi?= R x2sin 2 12sin 2 1 t PA 2121x (7.5) 116 ( ) ( ) ( ) ? ? ? ?? ? ??+?+?+? pi= R x2sin 2 12sin 2 1 t PB 2121x (7.6) ( ) ( ) ( ) ? ? ? ?? ? ????+ pi?= R y2cos 2 12cos 2 1r/rln2 t PA 2112y (7.7) ( ) ( ) ? ? ? ?? ? +??? pi?= R y2cos 2 12cos 2 1 t PB 21y (7.8) ??? ? ??? ? ?= ? yR xtan 1 1q (7.9) ??? ? ??? ? += ? yR xtan 1 2q (7.10) ( )221 yRxr ?+= (7.11) ( )222 yRxr ++= (7.12) P: load, t : thickness, and R : radius. Theoretical displacement fields calculated from equation (7.4) are illustrated in Figure 7.2 (b) and (c) in the form of contour maps, or moir? fringe patterns. The diameter and thickness of the disc were 25.4 mm and 3.17 mm, respectively. The applied load was 1800 N and the elastic constants for aluminum (E = 70 GPa, ? = 0.33) were used for the calculations. 117 (a) (b) (c) Figure 7.2 (a) Schematic diagram of a circular disc in compression, and (b) U and (c) V displacement fields obtained from a theoretical solution 118 7.1.3 Over-deterministic inverse approach 7.1.3.1 Displacement fields of circular disc in diametrical compression Equation (7.4) shows that E and ? are coupled non-linearly. Therefore, E and ? can be obtained simultaneously using either the U or the V displacement field. This is a very important advantage of the proposed approach. Numerous simple optical schemes are available for two-beam moir? interferometry (single field). Then which field should be used? Further investigation of equations (7.5) to (7.8) reveals that Ax and Bx are similar in magnitude but By is generally much smaller than Ay. Consequently, the V field displacements are much more sensitive to E than ?. The theoretical displacement fields of Figure 7.2 were utilized to investigate the sensitivity of displacement fields to elastic constants. The displacement fields were calculated for various combinations of E and ?; each variation was ?25%. The U displacements along the horizontal centerline (OA in Figure 7.2 b) and the corresponding V displacements along the vertical centerline (OB in Figure 7.2 c) are plotted in Figure 7.3 (a) and (b), respectively. The displacements were normalized by the displacement of the reference case at point A. The plots clearly indicate that the U displacements are uniquely defined for different combinations of E and ?, whereas the V displacements are virtually insensitive to ?. Only U displacements can be used for simultaneous measurement of both constants. In the following analyses, the term ?displacement? represents the U displacements for simplicity. 119 0.00 0.25 0.50 0.75 1.00 1.25 1.50 0.00 0.20 0.40 0.60 0.80 Normalized displacement along OA (x/R) Normalized Displacement ?,? ?,0.75? 0.75?,? 0.75?,1.25? 0.75?,0.75? 1.25?,? 1.25?,1.25? (a) 0.00 0.25 0.50 0.75 1.00 1.25 1.50 0.00 0.20 0.40 0.60 0.80 Normalized displacement along OB (y/R) Normalized Displacement ?,? ?,0.75? 0.75?,? 0.75?,1.25? 0.75?,0.75? 1.25?,? 1.25?,1.25? (b) Figure 7.3 Sensitivity of displacement fields to elastic constants, E and n. 120 7.1.3.2 Nonlinear least squares analysis The least-squares method has been used in a regression analysis [139]~[142]. The basic assumption that underlies this approach is that there are always differences between experimental results and theoretical values. Their relationship can be expressed as sru(x,y)U(x,y)e(x,y)e(x,y)=++ (7.13) where u and U are the theoretical and experimentally-measured displacements at each point (x,y), respectively, and es is the systematic displacement error and er is the random noise. Under an idealized condition where es and er are negligible, E and ? can be obtained using two arbitrary points in the U displacement field. In practice, the errors are not negligible and this is the rationale of the least-squares approach to fit the experimentally determined displacements to the theoretical solution. There are two types of random noise. One is random electrical noise introduced by the CCD camera and image grabber used for imaging and digital image processing. The other is defects or imperfections of the specimen grating. In practice, it is difficult to control accidental rigid-body motions of the specimen while applying external loads. Systematic errors are associated with these rigid-body motions. When studying deformations, we are interested in only relative displacements and the rigid-body displacements are inconsequential. For the proposed approach, however, absolute displacements are required and displacements altered by rigid-body motions must be considered. Rigid-body translation in the x direction changes the fringe order of every point in the U field by a constant, while rigid-body translation in the y or z direction does not affect the U displacements. Rigid-body in-plane rotation (about the z axis) alters the apparent 121 U displacements and it changes the fringe order by a linear function of y. Although out- of-plane rigid-body rotation about the x axis does not alter the fringe pattern, out-of-plane rigid-body rotation about the axis parallel to the grating lines (about the y axis) causes a foreshortening effect of the specimen grating [135]. This effect produces an apparent compressive strain of the specimen and thus apparent U displacement, which changes the fringe order by a linear function of x. Considering the above and using equation (7.1), the error function of least squares method can be expressed as ( ) ( ) 2M ii xx012r i1 s 2M iii xxx012r i1 s 1SuNxye 2f 11ABNxye EE2f = = ??=?++?+?+?+?? ?? ???=?+++?+?+?+?? ?? ? ? (7.14) where ixN is the fringe order at the point i in the fringe pattern; ixu is the corresponding theoretical displacement; M is the number of points used in the calculation; 0d is the fringe order change caused by the rigid body translation in the x direction; ( )1x? is the fringe order change caused by the out-of-plane rigid-body rotation about the y axis; and ( )2y? is the fringe order change caused by the in-plane rigid body translation about the z axis. 12and?? are fringe gradients, expressed as fringes per mm. Let E 1=? and E ??= (7.15) Then ( ) ( ) 2M iii xxx012r i1 s 1SABNxye 2f= ??=??+?++?+?+?+?? ??? (7.16) 122 The least squares criterion for the five independent variables requires ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? =?? =?? =?? =?? =?? 0S 0S 0S 0S 0S 2 1 0 d d d b a (7.17) This yields ( ) ( )( ) ( )( ) ( )( ) ( )( ) ( ) ( )( ) ( )( ) ( )( ) ( )( ) ( ) ( )( ) ( )( ) ( )( ) ( )( ) ( ) ( )( ) ( )( ) ( )( ) ( )( )??? ? ? ? ? ? ?? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ??? ? ??? ? ??? ? ??? ? ??? ? ??? ? = ?? ? ? ?? ? ? ? ?? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ??? ? ??? ? ??? ? ??? ? ??? ? ??? ??? ??? ? ??? ? ??? ? ??? ? ??? ? ??? ??? ??? ? ??? ? ??? ? ??? ? ??? ? ??? ??? ??? ??? ? ? ? ? ? ????? ????? ???? ????? ????? = = = = = ===== ===== ==== ===== ===== M 1i ii x 2 s M 1i ii x 2 s M 1i i x 2 s M 1i i x i x s M 1i i x i x s 2 1 0 M 1i 2i 2 s M 1i ii 2 s M 1i i 2 s M 1i ii x s M 1i ii x s M 1i ii 2 s M 1i 2i 2 s M 1i i 2 s M 1i ii x s M 1i ii x s M 1i i 2 s M 1i i 2 s 2 s M 1i i x s M 1i i x s M 1i ii x s M 1i ii x s M 1i i x s M 1i 2i x M 1i i x i x M 1i ii x s M 1i ii x s M 1i i x s M 1i i x i x M 1i 2i x yNf21 xNf21 Nf21 NBf21 NAf21 yf21yxf21yf21yBf21yAf21 yxf21xf21xf21xBf21xAf21 yf21xf21Mf21Bf21Af21 yBf21xBf21Bf21BBA yAf21xAf21Af21BAA (7.18) 123 The values of ?, ?, ?0, ?1 and ?2 can be determined by solving equation (7.18), and E and ? from equation (7.15). 7.1.4 Technical considerations for implementation Two technical issues arise when the proposed method is implemented: (1) the optimum number of data points for least-squares calculation and (2) the effect of the geometric center of the disc. A computer simulation was conducted to address these issues. 7.1.4.1 Systematic error and random noise The fringe patterns shown in Figure 7.2 were modified again to incorporate the systematic errors and the random noise. First, the maximum value of 0? was estimated to be 0.25 fringe order. This systematic error occurs when a zero fringe order is assigned to the center point of the fringe pattern, and 0 0.25fringeorder?= is a very conservative estimate. The maximum value 2? was chosen as 3.94 x 10-2 fringes per mm, which is equivalent to 1.0 fringe order change across the vertical diametric centerline (BB?). A moir? interferometer is usually equipped for adjustment of in-plane rotation, and the 1.0 fringe order difference across the diameter is readily discernable during adjustment. The value of 1d is directly related to the accidental out-of-plane rotation. In practice, the angle is determined by observing the spot of light reflected back to the source [135]. The magnitude of 1d can be determined by 124 2 1s df FL ???=? ?? ?? (7.19) where FL is the focal length of the collimating lens and d is the movement of the spot of light at the focal plane of the imaging lens. Considering a collimating lens with FL = 6? (150 mm) and d = 1 mm, 1? becomes 5.33 x 10-2 fringes per mm. The resultant U field fringe pattern with the above systematic noise is shown in Figure 7.4 (a). (a) (b) Figure 7.4 (a) Theoretical fringe patterns with systematic errors caused by rigid- body displacements, and (b) the corresponding fringe patterns after adding random noise with a maximum variance of 0.25 fringe. 125 The fringe pattern was then modified with random noise, which was specified as max r r s Neran() 2f= (7.20) where maxrN is the maximum variance of the random error in fringe order and ran() was random numbers ranging from ?1 to 1. The results of these modifications are illustrated in Figure 7.4 (b) for maxrN0.25= . The size of the digital fringe image is 1000 pixels by 1000 pixels, which is a typical resolution provided by a 1-inch format CCD camera. 7.1.4.2 Optimum number of data points It is desired to use all the points on the specimen for the nonlinear least-squares calculation because more data points can suppress random noise better and thus make the variables converge faster and more accurately. However, the theoretical solution is based on a point load, whereas the real experimental loading is applied over a small area (distributed load). The theoretical displacements near the loading area can be different from the experimental data because of the different loading condition and plastic deformations. In addition, the specimen grating along the free-edges can be damaged during replication, which can disturb the displacements along the free-edges. For these reasons, the points located in the central part of the specimen were used for least-squares calculation, as illustrated in Figure 7.5. The following simulation was conducted to determine the minimum number of data points that provides sufficient accuracy 126 Figure 7.5 Schematic diagram of the areas used in the nonlinear least squares calculation. The fringe pattern of Figure 7.4 (a) was first used to determine ? and ? with various values of m. The maximum value of m was limited to ?1? to avoid the points that are close to the loading areas and disc edges. All five variables converged within a few iterations and the values remained the same regardless of m. In order to evaluate the effect of random noise, five different sets of fringe patterns with maxrN0.25= were generated. The results with various values of m are plotted in Figure 7.6, where the values of ? and ? are normalized by those obtained from the fringe pattern of Figure 7.4 (a). The values of ? and ? converged within ?1% at m = 0.6. In typical moir? patterns, max rN does not exceed 0.25 fringe order and it is reasonable to assume that m = 0.6 will be sufficient for accurate determination of ? and ?, which corresponds to 90,000 data points when the digital fringe image has a resolution of 1000 pixels by 1000 pixels. 127 0.80 0.90 1.00 1.10 1.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 m Normalized a Set 1 Set 2 Set 3 Set 4 Set 5 (a) 0.80 0.90 1.00 1.10 1.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00 m Normalized b Set 1 Set 2 Set 3 Set 4 Set 5 (b) Figure 7.6 Effect of calculation area on a and b. 128 7.1.4.3 Effect of geometric center of the disc It is difficult to locate the exact location of the origin in the digital fringe patterns. This also produces a systematic error in the experimentally determined displacement fields. Mathematically, the error, cpe , can be expressed as ( ) ( )cp00eUx,yUxx,yy=??? (7.21) where 00xandy are the x and y coordinates of an incorrect origin. It is important to note that this error is different from that caused by the horizontal rigid-body translation ( 0? ). The 0? changes the fringe order uniformly and thus the relative displacement between any two points in the field remains unchanged. However, the change in fringe order produced by cpe varies from point to point and the relative displacement changes with 00xandy . A computer simulation was conducted to illustrate the effect of cpe . In the simulation, the position of origin was intentionally moved to a point of (0.02R, 0.02R) and the values of ? and ? were determined with m greater than 0.6. The results are shown in Figure 7.7. The values of ? and ? were underestimated significantly and converged slowly as m increased. Another simulation with an origin of (?0.02R, ?0.02R) was also conducted and the results were virtually identical to the plot shown in Figure 7.7. The simulation indicates that only the correct location of the origin can provide stable outputs as m increase from 0.6 to 1.0. 129 0.80 0.90 1.00 1.10 1.20 0.60 0.70 0.80 0.90 1.00 m Normalized a Set 1 Set 2 Set 3 Set 4 Set 5 (a) 0.80 0.90 1.00 1.10 1.20 0.60 0.70 0.80 0.90 1.00 m Normalized b Set 1 Set 2 Set 3 Set 4 Set 5 (b) Figure 7.7 Effect of location of the origin on a and b. 130 The 00xandy can be considered as variables in the nonlinear least-squares algorithm. Unlike other systematic errors, however, they are coupled with the material constants non-linearly. This would affect the stable convergence of the constants. Instead, a simple solution using a concept of minimum variance is proposed to find the correct location of the origin. Recalling the analysis shown in the previous section, the values of ? and ? converge and remain unchanged if a region with m > 0.6 is used for calculation in spite of the random noise. In other words, the variance in ? and ? for m greater 0.6 should provide minimum variance if the correct origin is used in the calculation. This approach can be expressed mathematically as ? ? ? ? ? ? ? ? ??? ? ??? ? ?+? ? ?? ? ? ?= ?? == 2K 1i i 2K 1i i ? ? ? ? 2 1Var b bb a aa (7.22) where K i i1 K i i1 1? K 1? K = = ?=? ?=? ? ? (7.23) and K is the number of areas with m > 0.6 used for the calculation, ia and ib are the values calculated from equation (7.18) for the i-th calculation area. A point with the minimum variances of ? and ? is selected as the origin, and the E and ? are determined from: 1E ? ? ? = ? ??= ? (7.24) 131 7.1.5 Experimental verification A moir? experiment was conducted to verify the proposed method. The specimen was a circular disc, shown schematically in Figure 7.8 (a); it was made of 6061-T6 aluminum alloy. The diameter was 25.4 mm and the thickness was 3.17 mm. Two spherical indents with a radius of 3.2 mm were made and the load was transferred through steel balls to ensure accurate alignment of diametrical compression. A general- purpose compression loading fixture was utilized for the experiment. As shown schematically in Figure 7.8 (b), the compressive displacements were achieved by driving the lower wedge to the left. The load was measured with an accuracy of 0.01 N by an electrical load cell. Figure 7.8 (c) and (d) show the photos of the real experiment setup. The U field fringe patterns at 1800 N are shown in Figure 7.9; (a)-(d) is a series of four phase-shifted fringe patterns for this load level; (e) and (f) are the corresponding wrapped and unwrapped phase maps. The fringe patterns at a higher load level (2200 N) are shown in Figure 7.10. The proposed algorithm was implemented with these fringe patterns. First, the geometric center of the image was marked as an initial location of the origin. Then the area of possible locations of the true origin was selected as 0.04R by 0.04R (n = 0.04 in Fig. 4), which corresponded to an area of 20 pixels by 20 pixels in the digital fringe patterns. The values of ? and ? were calculated for each point within the area by equation (7.18) using m > 0.6 with an interval of m = 0.1. The variances were calculated by equation (7.22) and the point that produced the minimum variance was chosen as the true origin. 132 The results obtained using the correct origin are plotted in Figure 7.11 for m > 0.6, where ia and ib were normalized by ?? and??, respectively. The material constants were evaluated by equation (7.24) and the results are summarized in Table 7.1. The results agree well with the values of the handbook [143], which confirms the validity of the proposed method. (a) (b) 133 (c) (d) Figure 7.8 (a) Schematic diagram of the specimen; (b) Schematic diagram of compression loading fixture used in the moir? experiments; (c) Photo of the experiment setup and (d) the specimen loaded in the loading fixture. 134 (a) 0 (b) pi/2 (c) pi (d) 3pi/2 (e) wrapped phase map (f) unwrapped phase map Figure 7.9 U field fringe patterns obtained at 1800 N; (a)-(d) four phase-shifted moir? patterns; (e) and (f) phase maps before and after unwrapping process. 135 (a) 0 (b) pi/2 (c) pi (d) 3pi/2 (e) wrapped phase map (f) unwrapped phase map Figure 7.10 U field fringe patterns obtained at 2200 N; (a)-(d) four phase-shifted moir? patterns; (e) and (f) phase maps before and after unwrapping process. 136 0.80 0.90 1.00 1.10 1.20 0.60 0.70 0.80 0.90 1.00 m Normalzied a 1800 N 2200 N (a) 0.80 0.90 1.00 1.10 1.20 0.60 0.70 0.80 0.90 1.00 m Normalized b 1800 N 2200 N (b) Figure 7.11 Values of a and b obtained from the fringe patterns in Fig. 7.8 and 7.9. 137 Table 7.1 Experimental results E (GPa) ? 1800 N 70.4?1.1 0.32?0.01 2200 N 69.7?1.1 0.31?0.01 Handbook [143] 69.7 0.33 7.1.6 Discussion The proposed method utilizes only the U field fringe pattern to determine the two elastic constants. A variety of simple optical configurations of two-beam moir? interferometry can be employed. Rigid-body displacements induced by practical loading systems cannot be avoided. To account for these effects, the magnitudes of 1? and 2? were 3.94 x 10-2 and 5.33 x 10-2 fringes per mm, respectively. These displacement gradients are equivalent to only 17 and 22 micro strain, which are negligible for most deformation analyses. Together with the unavoidable random noise, however, it was demonstrated that these small displacement gradients could produce significant uncertainties in the elastic constants. With the high accuracy provided by moir? interferometry and with the over-deterministic approach, the proposed algorithm handles the rigid-body motions effectively. The use of simple optical configurations to record a single displacement field, together with the powerful algorithm that is capable of handling the accidental rigid-body displacements, provides a practical method for advanced engineering materials. 138 7.1.7 Summary An inverse method has been proposed to determine the elastic constants E and n. The method is based on the theoretical displacement of a circular disc in diametrical compression. High sensitivity moir? interferometry combined with the phase-shifting technique provides the experimental data in the form of a full-field displacement map. An over-deterministic approach using the nonlinear least squares method is implemented to fit the experimentally determined displacements to the theoretical solution. A computer simulation was conducted to investigate the effect of systematic errors and random noise encountered in implementing the proposed method. The optimum number of data points for least-squares calculation was determined as 0.36R2 of the number of pixels in the image. It was found that the results were sensitive to the location of the geometric center, and a simple procedure using minimum variance was utilized to identify the correct origin. A moir? experiment was conducted to verify the proposed method experimentally. With only one moir? field, two elastic constants, E and n, were determined simultaneously with high accuracy. 139 7.2 In-plane deformation and strain measurement of electronic packaging components Moir? interferometry is the most widely used experimental technique for thermal deformation analyses of electronic packaging [144]~[153]. The typical sensitivity of moir? interferometry is 0.417?m. As the components and structures are made smaller and smaller, microscopic moir? interferometry, which can provide higher resolution and sensitivity, has also become important. 7.2.1 In-plane deformation measurement using moir? interferometry Moir? interferometry uses a high frequency diffraction grating. Because of the defect and imperfection of the high-frequency specimen grating, the fringe images obtained from moir? interferometry usually contain much more noise than those obtained in out- of-plane displacement measurement methods. Figure 7.12 shows the experiment results of the deformations of a plastic quad flat package (PQFP) due the thermal loading and hygroscopic swelling [154]. PQFP is a plastic encapsulated microcircuit which consists of a silicon chip, a metal support or leadframe, wires than electrically attach the chip?s circuits to the leadframe, and a plastic epoxy encapsulating material, or mold compound, to protect the chip and the wire interconnects. The mold compound is a composite material made of an epoxy matrix that encompasses silica fillers, stress relief agents, flame-retardants, and many other additives. In spite of many advantages, one important disadvantage of PEMs is that the polymeric mold compound absorbs moisture when exposed to a humid environment due to the 140 polymer-water affinity action. The goal of the experiment in Figure 7.12 was to investigate the effect of hygroscopic swelling on the package and to compare it with that of thermal expansion. (a) Fringe patterns generated due to thermal loading ?T = -60?C (b) Fringe patterns generated due to hygroscopic swelling under a virtual saturation state at %85RH Figure 7.12 Moir? fringe pattern of PQFP package with hygroscopic and thermal deformations Figure 7.13 shows the deformed configuration. It is obvious that hygroscopic swelling brings much more deformation than a thermal loading (?T = -60?C), thus 141 hygroscopic swelling must be considered for accurate reliability assessment when PQFP packages are subjected to environments where the relative humidity fluctuates. It is also worth noting that the deformed shape due to hygroscopic swelling is convex whereas the deformed shape due to heating is concave (convex for cooling). (a) Deformed shape due to hygroscopic swelling (b) Deformed shape due to thermal loading (cooling) Figure 7.13 2-D deformations (Deformations are magnified with same magnification) 142 7.2.2 In-plane deformation and strain measurement using microscopic moir? interferometry Microelectronics devices contain many very tiny electronic components within an active silicon chip, such as transistors, capacitors, resistors, etc. A silicon chip requires protection from the environment as well as electrical and mechanical connections to the surrounding components. The various conducting and insulating materials involved in the devices have different coefficients of thermal expansions (CTEs). When electrical power is applied, the device is subjected to a temperature excursion and each material expands at a different rate. This non-uniform CTE produces thermally induced mechanical stresses within the device assembly. As the components and structures involved in high-end microelectronics devices are made smaller, the thermal gradient increases and the strain concentrations become more serious. Hence, there is a continuously increasing activity in experimental analysis, both for specific studies and for guidance of numerical analyses. Moir? interferometry and microscopic moir? interferometry are the leading methods for experimental analyses [2]~ [4][121][155]. One of several purposes of a chip carrier is to provide conducting paths between the extremely compact circuits on the chip and the more widely spaced terminals on the PCB. The micro via technology enabled the industry to produce laminate substrates with a high density, and fine pitch conductors, called surface laminar circuit (SLC), as required for advanced assemblies. Extensive research and development efforts have been and are being made to perfect the underfill process for the high-density organic substrates, and to develop optimum underfill materials for the larger silicon devices. An 143 important trend in newly developed underfill materials is its increased Young?s modulus, which increases the coupling between the silicon chip and the substrate. This high degree of coupling transfers the CTE mismatch-induced-loading to the build-up layers of the substrate. Figure 7.14 (a) shows the schematic diagram of the flip-chip assembly on a high-density substrate used in reference [156] to quantify the thermal-mechanical deformations of the microstructures within the build-up layers. Microscopic moir? interferometry [157] was employed in the experiment. Two specimen configurations were analyzed to study the deformations induced by the subsequent manufacturing process: a bare substrate and a flip-chip package. In the assembly, a silicon chip was attached to a high-density substrate by solder bumps and the gap between the solder bumps was filled with an underfill [156]. Figure 7.14 Schematic diagram of the flip-chip assembly on a high-density substrate (a) before and (b) after specimen preparation 144 The specimens were cut and ground to expose the desired microstructures as illustrated schematically in Figure 7.14 (b), where the inset depicts the detailed microstructures within the build-up layer. The specimen grating was replicated at 92?C and the fringes were recorded at room temperature of 22?C. Therefore, the thermal loading is ?T = -70?C. The displacement fields for a small region containing the microstructures were recorded by microscopic moir? interferometry. The region is marked by a dashed box in Figure 7.15; it is approximately 500 ?m by 375 ?m. The resultant fringe patterns are shown in Figure 7.16 (a) for the bare substrate and in Figure 7.16 (b) for the flip-chip assembly. O/DFM was employed to produce a displacement contour interval of 52 nm/fringe. Figure 7.15 Micrographs of the region of interest. 145 Figure 7.16 Microscopic U and V displacement fields of (a) the bare substrate and (b) the flip-chip assembly. To investigate the effect of the chip and underfill, the deformations of the solder and metal via were analyzed using hybrid semi-automatic fringe processing technique. The deformed shape of the solder and the metal via in the substrate and the flip-chip assembly are plotted in Figure 7.17 and Figure 7.18, respectively. The deformed shapes were exaggerated with same magnification in both figures. It is evident that deformations of the solder and the metal via increased significantly after the chip was assembled to the substrate. 146 (a) Undeformed shape (b) Deformed shape Figure 7.17 Deformations of the solder and the metal via in substrate. (a) Undeformed shape (b) Deformed shape Figure 7.18 Deformations of the solder and the metal via in Flip-chip assembly. The normal strain ?y strain distributions are plotted in Figure 7.19. Compared with the strain field in substrate, strain concentration happens in flip-chip assembly. The maximum strain in assembly is about three times of the maximum strain the substrate. 147 Figure 7.19 Strain ey distributions in the solder and the metal via. 7.2.3 Coefficient of thermal expansion (CTE) measurement using moir? interferometry Because coefficient of thermal expansion (CTE) mismatch can reduce the reliability of electronic packaging systems by causing localized stress, CTE measurement is important for optimizing the packaging design. Real-time moir? interferometry has been successfully used to measure the CTE of the electronic packaging components [158]. Mathematically, CTE can be calculated by T? ??=? (7.25) where ? is CTE, ?? is the strain change due to a temperature change of ?T. Apparently, CTE measurement is actually a thermal strain measurement. Figure 7.20 and Figure 7.21 show some fringe patterns for the CTE measurement of the printed circuit board (PCB) and the module in a PBGA package. Compared with manual method which selects several lines across the specimen for the CTE calculation, 148 computer-aided fringe analysis can calculate the whole-field CTEs and the average CTE can be regarded as the specimen?s CTE. Table 7.2 shows the results of CTEs of PCB and module. (a) U at 60 ?C (b) V at 60 ?C (c) U at 100 ?C (d) V at 100 ?C Figure 7.20 CTE measurement of PCB. 149 (a) U at 60 ?C (b) V at 60 ?C (c) U at 100 ?C (d) V at 100 ?C Figure 7.21 CTE measurement of module. Table 7.2 Results of CTE measurement at 60~100 C CTE: X direction (ppm/?C) CTE: Y direction (ppm/?C) PCB 18.0 17.5 Module 6.8 6.7 150 8. CONCLUSIONS 8.1 Conclusions Photomechanics methods produce fringe patterns, which generally provide the full field information of displacement field. Traditionally, the fringe patterns have been analyzed manually for experiment analyses. As digital image processing has become more accessible, computer-aided automatic fringe pattern analyses have become practical. Although numerous image-processing algorithms have been developed to complement interferometric measurement techniques, their extensions to general fringe analyses have been limited because of inherent optical noise encountered in the modern photomechanics techniques. To enhance the applications of photomechanics methods to advanced engineering problems, a robust general-purpose computer-aided fringe analysis tool was developed in this dissertation. The major contributions made in this dissertation include: (1) Complete survey and investigation on the existing fringe image processing schemes; the most appropriate image processing schemes and their limitations were identified. (2) Development of self-adaptive fringe filtering algorithm. This algorithm is a spatial low-pass filtering. Unlike the conventional digital image filtering, the filter window in this algorithm is selected in a self-adaptive manner; it considers not only the orientations of the local fringes but also the local fringe densities. 151 This algorithm is very effective for suppressing and reducing noise without blurring and losing useful fringe information. (3) Quantitative evaluation of the automatic fringe analysis techniques: Fourier transform technique and phase shifting technique. Although the techniques are widely used, the evaluation shows that neither of the two techniques can handle discontinuities effectively. The investigation also shows that these two techniques are effective for determining fractional fringe orders, but should not be used for calculation of fringe gradient. (4) Development of enhanced random phase shifting algorithm. The conventional phase shifting algorithm requires pre-determined or known phase shift amounts. The random phase shifting algorithm is a least squares iteration procedure; it can detect the phase shift amounts and the full-field phase distributions automatically and simultaneously. Using the random phase shifting algorithm, the phase shift amounts can be any random values; thus an accurate phase shifting devices are not required in the experiment. (5) Development of semi-automatic fringe order assignment and interpolation. A series of algorithms for fringe order assignment and interpolation are proposed and implemented, which include semi-automatic fringe order assignment, 1-D segment-by-segment curve fitting interpolation, improved 1-D fringe order interpolation considering fringe orientation, and 2-D tile-by-tile interpolation. These algorithms make the existing fringe centering technique practical and increase applicability of the fringe centering technique significantly. 152 (6) Development of hybrid O/DFM (Optical/digital fringe multiplication) fringe centering technique. This hybrid technique combines the advantages of the O/DFM, the phase shifting technique and the fringe centering technique. The hybrid O/DFM fringe centering technique is a very effective semi-automatic processing technique and it can be employed to obtain full-field fractional fringe orders (displacements) and their gradients (strains) accurately. (7) Development of a Windows GUI-based expert system (software) for interferogram fringe analysis and processing. The software system includes the general digital image processing algorithms and all the algorithms introduced in this dissertation. (8) Application using infrared diffraction interferometer for the co-planarity measurement of high-density solder bump patterns. The method utilizes long wavelength (? ?= 10.6 ?m), coherent infrared laser light, which serves to reduce the apparent roughness of test objects, and enables the regularly spaced solder bump arrays to produce well-defined diffracted wavefronts. The expert system was utilized to produce a full-field surface topography map. (9) Development of an inverse method to determine elastic constants using circular disc and moir? interferometry. This inverse method uses a circular disc in diametrical compression to simultaneously determine the two elastic constants, E and n, from a single displacement map. Moir? interferometry combined with the phase-shifting technique provides a full-field displacement field. An over- deterministic approach using the nonlinear least squares method is implemented to fit the experimentally determined displacements to the theoretical solution. 153 (10) Applications of photomechanics and computer-aided digital fringe processing to electronic packaging. The computer-aided fringe processing techniques have been successfully applied to obtain the out-of-plane shape and warpage, in-plane deformation and strain of various electronic and microelectronic packaging components. 8.2 Future work This dissertation work presents a variety of robust computer-aided fringe analysis techniques for the general?purpose photomechanics fringe analysis. The expert system can be used to obtain whole field fringe orders (usually, displacements) with high accuracy automatically; however, the gradient calculation is very sensitive to noise and very small change of fringe orders can result in large gradient variation. Therefore, the whole-field gradient calculation is usually a qualitative analysis. Fortunately, the engineering problems usually require a quantitative analysis at one single point, along one single line or in one small area. With user?s judgment, such a quantitative analysis is possible in the expert system; nevertheless, more work is required to achieve quantitative analysis of fringe order gradients automatically. A least squares inverse approach was developed in this dissertation to determine the elastic constants using whole-field displacement information; a similar approach can be extended to measure the residual stresses using moir? hole drilling method. Moir? interferometry combined with computer-aided fringe analysis can provide whole-field displacement information around the moir? hole, an over-deterministic least squares 154 approach can then be employed to obtain the orientations and magnitudes of the residual stresses. This future work will provide a new way for residual stress measurement with high accuracy. 155 REFERENCES [1] B. Han, ?Recent advancement of moir? and microscopic moir? interferometry for thermal deformation analyses of microelectronics devices,? Experimental Mechanics, Vol. 38, No. 4, pp. 278?288, 1998. [2] B. Han and Y. Guo, ?Thermal deformation analysis of various electronic packaging products by moir? and microscopic moir? interferometry,? Journal of Electronic Packaging, Trans. ASME, Vol. 117, pp. 185?191, 1995. [3] B. Han, Y. Guo, C. K. Lim, and D. Caletka, ?Verification of numerical models used in microelectronics packaging design by interferometric displacement measurement methods,? Journal of Electronic Packaging, Trans. ASME, vol. 118, pp. 157?163, 1996. [4] P. Kunthong, K. Verma, and B. Han, ?Effect of underfill on C4 bumps and surface laminar circuit: an experimental study,? Proceedings of the 1999 IMAPS Conference, Chicago, IL, October 1999. [5] K. Verma, B. Han, S.-B. Park, and W. Ackerman, ?On the design parameters of flip-chip PBGA package assembly for optimum solder ball reliability,? IEEE Trans. Components and Packaging Technologies, Vol. 24, No. 2, pp. 300?307, 2001. [6] B. Han, ?Optical measurement of flip-chip package warpage and its effect on thermal interfaces,? Electronics Cooling, Vol. 9, No. 1, pp. 10?16, February 2003. 156 [7] K. Crennel, ?Introductory digital image processing?, In: Robinson D. W., Reid G. T., eds., Interferogram analysis. London: IOP Publishing, pp. 1-22, 1993. [8] T. Kreis, Holographic Interferometry. Berlin: Akademie Verlag, 1996 [9] D. Malacara, M. Servin and Z. Malacara, Interferogram analysis for optical testing. New York: Marcel Dekker Inc. 1998 [10] A. Jain, Fundamentals of digital image processing. London: Prentice-Hall International Inc. 1989 [11] G. Jin, Computer-aided optical measurements. Beijing, Tsinghua University Press. 1997 [12] T. Yatagai, ?Internsity based analysis methods,? In: Robinson D W, G T Reid, eds., Interferogram analysis. London: IOP Publishing, pp. 73~93, 1993. [13] G. Mastin and D. Ghiglia, ?Digital extraction of interference fringe contours,? Applied Optics, Vol. 24, pp. 1727-1728. 1985. [14] T. Chen and C. Taylar, ?Computerized fringe analysis in photomechanics?, Experimental Mechanics, Vol. 29, No. 3, pp. 323-329, 1989. [15] Y. Morimoto, ?Digital image processing,? In: Kobayashi A. eds., Handbook on experimental mechanics. 2nd edition. SEM, 1994. [16] T. Yatagai, S. Nakadate, M. Idesawa, et al. ?Automatic fringe analysis using digital image processing techniques,? Optical Engineering, Vol. 21, No. 3, pp. 432~435, 1982. 157 [17] T. Yatagai, ?Automatic fringe analysis techniques in Japan,? Optics and Lasers in Engineering, Vol. 15, No. 2, pp. 79~91, 1991. [18] A. Bastawros and A. Voloshin, ?Thermal strain measurements in electronic packages through fractional fringe moir? interferometry?, Journal of Electronic Packaging, Trans. ASME. Vol. 112, No. 4, pp. 303-308, 1990. [19] F. Bertolino and F. Ginesu, ?Semiautomatic approach of grating techniques,? Optics and Lasers in Engineering, Vol. 19, pp. 313-323, 1993. [20] A. Ennos, D. Robinson and D. Williams, ?Automatic fringe analysis in holographic interferometry,? Optica Acta, Vol. 32, pp. 135-145, 1985. [21] S. Krishnaswamy, ?Algorithm for computer tracing of interference fringes?, Applied Optics, Vol. 30, pp. 1624-1628, 1991. [22] J. Aparicio, ?Improved algorithm for the analysis of holographic interferograms,? Optical engineering, Vol. 32, No. 5, pp. 963-969, 1993. [23] D. Zhang, M. Ma and D. Arola, ?Fringe skeletonizing using an improved derivative sign binary method,? Optics and Lasers in Engineering, Vol. 37, pp. 51-62, 2002. [24] D. Zhang, D. Arola and M. Ma, ?Digital image processing of moir? fringe patterns with an application to fractures in bovine dentin?, Optical engineering, Vol. 41, No. 5, pp. 1115-1121, 2002. [25] B. Han, ?Interferometric methods with enhanced sensitivity by optical/digital fringe multiplication,? Applied Optics, Vol. 32, pp. 4713-4718, 1993. 158 [26] X. He, J. Dai and S. Liu, "Automatic warpage measurement of advanced packages by optical method," ASME, EEP-Vol.24, Thermo-Mechanical Characterization of Evolving Packaging Materials and Structures, Nov., 15-20 , Anaheim, California, pp. 65-70, 1998. [27] G. Johnson, D. Leiner and D. Moore, ?Phase-locked interferometry,? Optical Engineering, Vol. 18, pp. 46-52, 1979. [28] A. Rosenfeld, ?A characterization of parallel thinning algorithms?, Information and Control, Vol. 29, pp. 286-291, 1975. [29] A. Rosenfeld and A. Kak, Digital Picture Processing, Vol. I, 2nd edition, Academic Press, New York, 1982 [30] C. Hilditch, ?Comparison of thinning algorithms on a parallel processor?, Image and Vision Computing, Vol. 1, 115-132, 1983. [31] X. Wang, C program for image processing. Hefei: CSTU press, 1994 [32] U. Mieth and W. Osten, ?Three methods for the interpolation of phase values between fringe pattern skeletons?. In Osten W, Pryputniewicz R., , Reid G., Rottenkolber H. eds., Fringe '89, Automatic Processing of Fringe Patterns. Physical research, Akademie-Verlag Berlin, pp. 118-123, 1989. [33] J. Schemm and C. Vest, ?Fringe pattern recognition and interpolation using nonlinear regression analysis,? Applied Optics, Vol. 22, pp. 2850-2853, 1983. [34] R. Muller and L. Saackle, ?Complete automatic analysis of photoelastic fringes,? Experimental Mechanics, Vol. 19, No. 7, pp. 245~251, 1979. 159 [35] A. Colin and W. Osten, ?Automatic support for consistent labeling of skeletonized fringe patterns,? J. Mod. Opt. Vol. 42, pp. 945-954, 1995. [36] D. Burton and M. Lalor, ?Managing some of the problems of Fourier fringe analysis,?. In Reid G., ed., Fringe Pattern Analysis, Proc. of Soc. Photo-Opt, Instr. Enq., 1163, pp. 149-160. 1989. [37] W. Macy, ?Two-dimensional fringe pattern analysis,? Applied Optics, Vol. 22, pp. 3898-3901, 1983. [38] P. Ransom and J. Kokal, ?Interferogram, analysis by a modify sinusoid fitting technique,? Applied Optics, Vol. 25, pp. 4199-4204, 1986. [39] Q. Yu, X. Liu and X. Sun, ?Double-image and single-image phase shifting methods for phase measurement,? Optik, Vol. 109, No. 2, pp. 89~93, 1998. [40] K. Freischlad and C. Koliopoloulos, ?Fourier description of digital phase- measuring interferometry,? Journal of the Optical Society of America A. Vol. 7, pp. 542-551, 1990. [41] T. Kreis, ?Digital holographic interference-phase measurement using the Fourier transform method,? Journal of the Optical Society of America A. Vol. 3, pp. 847- 855, 1986. [42] T. Kreis, ?Fourier-transform evaluation of holographic interference patterns,? In F.-P. Chiang, ed., International Conference on Photomechanics and Speckle Metroloqy. Proc. of Soc. Photo-Opt. Instr. Eng., 814, pp. 365-371, 1987. 160 [43] A. Malcolm and D. Burton, ?The relationship between Fourier fringe analysis and the FFT,? In Prypntniewicz R., ed.. Laser Interferometry IV: Computer-Aided Interferometry. Proc. of Soc. Photo-Opt. Instr. Eng. 1553. pp. 286-297, 1991. [44] Q. Ru, T. Honda, J. Tsujinchi and N. Oliyama., ?Fringe analysis by using 2-D- Fresnel transforms,? Optics Communications, Vol. 66, pp. 21-24, 1988. [45] M. Takeda, H. Ina and S. Kobayashi, ? Fourier-transform method of fringe- pattern analysis for computer-based topography and interferometry,? Journal of the Optical Society of America A., Vol. 72, pp. 156-160, 1982. [46] T. Kreis, ?Automatic evaluation of interference patterns.? In W. Juptner, ed., Holography Techniques and Applications. Proc. of Soc. Photo-Opt. Instr. Ena., 1026. pp. 80-89, 1988. [47] P. Zanetta, D. Albrecht, G. Schirripa and D. Paoletti D, ?Application of fast fourier transform techniques to the quantitative analysis of holographic and TV- holographic interferograms,? Optik, Vol. 97, pp. 47-52. 1994. [48] D. Bone, H. Bachor and R. Sandeman, ?Fringe-pattern analysis using 2-D Fourier transform,? Applied Optics, Vol. 25, pp. 1653- 1660, 1986. [49] P. Bryanston-Cross, C. Quan and T. Judge, ?Application of the FFT method for the quantitative extraction of information from high-resolution interferometric and photoelastic data,? Optics and Laser Technology, Vol. 26, pp. 147-155, 1994. [50] C. Gorecki, ?Interferogram analysis using a Fourier transform method for automatic 3D surface measurement.? Pure Appl. Opt., Vol. 1, pp. 103-110, 1992. 161 [51] J. Gu and F. Chen, ?Fast Fourier transform, iteration, and least-squares-fit demodulation image processing for analysis of single-carrier fringe pattern?. Journal of the Optical Society of America A, Vol. 12, pp. 2159-2164, 1995. [52] K. Nugent, ?Interferogram analysis using an accurate fully automatic algorithm?. Applied Optics, Vol. 24, pp. 3101-3105, 1985. [53] C. Roddier and F. Roddier, ?Interferogram analysis using Fourier transform techniques.? Applied Optics, Vol. 26, pp. 1668-1673, 1987. [54] D. Tentori and D. Salazar, ?Hohogram interferometry: carrier fringes,? Applied Optics. Vol. 30, pp. 5157-5158, 1991. [55] K. Womack, ?Frequency domain description of interferogram analysis,? Optical Engineering, Vol. 23, pp. 369-400, 1984. [56] R. Green, J. Walker and D. Robinson, ?Investigation of the Fourier-transform method of fringe pattern analysis,? Optics and Lasers in Engineering, Vol. 8, pp. 29-14, 1988. [57] B. Breuckmann and W. Thieme, ?Computer-aided analysis of holographic interferograms using the phase-shift method,? Applied Optics, Vol. 24, pp. 2145- 2149, 1985. [58] G. Stoilov, T. Dragostinov, ?Phase-stepping interferometry: five-frame algorithm with an arbitrary step,? Optics and Lasers in Engineering, Vol. 28, No. 1, pp. 61~69, 1997. 162 [59] D. Robinson, D. Williams, ?Digital phase stepping speckle interferometry,? Optics Communications, Vol. 57, No. 1, pp. 26~30, 1986. [60] G. Kaufmann, ?Phase shifting in speckle photography,? Optics and Lasers in Engineering, Vol. 29, No. 2/3, pp. 117~124, 1998. [61] K. Creath, ?Temporal phase measurement methods,? In: Robinson D W, G T Reid, eds. Interferogram analysis. London: IOP Publishing, pp. 95~140, 1993. [62] M. Chang, C. Hu, P. Lam and J. Wyant, ?HigH precision deformation measurement by digital phase shifting interferometry,? Applied Optics, Vol. 24, pp. 3780-3783, 1985. [63] Y. Cheng and J. Wyant, ?Phase shifter calibration in phase-shifting interferometry,? Applied Optics, Vol. 24, pp. 3049-3052, 1985. [64] K. Kinnstaetter, A. Lohmann, J. Schwider and N. Streibl, ?Accuracy of phase shifting interferometry,? Applied Optics, Vol. 27, pp. 5082-5089, 1988. [65] P. Groot, ?Derivation of algorithms for phase-shifting interferometry using the concept of a data-sampling window?. Applied Optics, Vol. 34, pp. 4723-4730, 1995. [66] J. Huntley, ?Automated fringe pattern analysis in experimental mechanics: a review?. Journal of Strain Analysis. Vol. 33, No. 2, pp. 105-125, 1998. [67] C. Poon, M. Kujawinska M and C. Ruiz, ?Automated fringe pattern analysis for moir? interferometry?, Experimental Mechanics, Vol. 3, pp. 234-241, 1993. 163 [68] C. Poon, M. Kujawinska and C. Ruiz, ?Spatial-carrier phase shifting method of fringe analysis for moir? interferometry?, Journal of Strain Analysis, Vol. 28, No. 2, pp. 79-88, 1993. [69] K. Perry and J. McKelvie, ?A comparison of phase shifting and Fourier methods in the analysis of discontinuous fringe patterns,? Optics and Lasers in Engineering, Vol. 19, pp. 269-284, 1993. [70] G. Lai and T. Yatagaib, ?Generalized phase-shifting interferometry,? Journal of the Optical Society of America A. Vol. 8, pp. 822-827, 1991. [71] T. Kreis, J. Geldmacher and W. Juptner, ?A comparison of interference phase determination methods with respect to achievable accuracy,?. In Juptner W. and Osten W. eds., Fringe '93. pp. 51-59, 1993. [72] K. Okada, A. Sato and J. Tsujiuchi, ?Simultaneous calculation of phase distribution and scanning phase shift in phase shifting interferometry?, Optics Communication, Vol. 84, No. 3&4, pp. 118-124, 1991. [73] K. Perry and J. McKelvie, ?Refernce phase shift determination in phase shifting interferometry,? Optics and Lasers in Engineering, Vol. 22, pp. 77-90, 1995. [74] J. Greivenkamp, ?Generalized data reduction for heterodyne interferometry?. Optical Engineering, Vol. 23, pp. 350-352, 1984. [75] E. Vikhagen and O. Lokberg, ?Detection of defects in composite materials by television holography and image processing,? Materials Evaluation, Vol. 48, No. 2, pp. 244~248, 1990. 164 [76] J. Gu, ?Phase iteration and spacing iteration algorithms for speckle fringe pattern processing,? Optics and Lasers in Engineering, Vol. 29, No. 2/3, pp. 145~158, 1998. [77] V. Srinivasan, S. Yeo and Chaturvedi. ?Fringe processing and analysis with a neural network?. Optical. Engineering, Vol. 33, pp. 1166-1171, 1994. [78] D. Ghiglia and M. Pritt, Two-dimensional phase unwrapping: Theory, algorithms and software, John Wiley & Sons, Inc., New York 1998. [79] A. Spik and D. Robinson, ?Investigation of the cellular automata method for phase unwrapping and its implementation on an array processor,? Optics and Lasers in Engineering, Vol. 14, No. 1, pp. 25~37, 1991. [80] K. Creath, ?Phase measurement interferometry techniques,? Progress in Optics, XXIV, pp. 349~393, 1988. [81] M. Owner-Peterson, ?Phase-map unwrapping: a comparison of some traditional methods and a presentation of a new approach,? In: Proceeding of SPIE---The International Society for Optical Engineering Industrial Application of Holographic and Speckle Measuring Techniques. Bellingham: ISOE, pp. 73~82, 1991. [82] D. Robinson, ?Phase unwrapping methods,? In: Robinson D W, G T Reid, eds. Interferogram analysis. London: IOP Publishing, pp. 194~229, 1993. [83] H. Chang, C. Chen, C. Lee, et al. ?The tapestry cellular automata phase unwrapping algorithm for interferogram analysis,? Optics and Lasers in Engineering, Vol. 30, No. 6, pp. 487~502, 1998. 165 [84] J. Wingerden, H. Frankena, and C. Smorenburg, ?Linear approximation for measurement errors in phase shifting interferometry,? Applied Optics, Vol. 30, pp. 2718-2729, 1991. [85] D. Birton, and M. Lalor, ?Multichannel Fourier fringe analysis as an aid to automatic phase unwrapping.? Applied. Optics, Vol. 33, pp. 2939-2918, 1991. [86] D. Bone, ?Fourier fringe analysis: the two-dimensional phase unwrapping problem?. Applied Optics, Vol. 30, pp. 3627-3632, 1991. [87] K. Itoh, ?Analysis of the phase unwrapping algorithm,? Applied. Optics, Vol. 21, pp. 2470, 1982. [88] A. Baidi, ?Two-dimensional phase unwrapping by quad-tree decomposition?, Applied Optics, Vol. 40, No. 8, pp. 1187-1194, 2001. [89] P. Bryanston-Cross and C. Quan, ?Examples of automatic phase unwrapping applied to interferometric and photoelastic images?. In Juptner W. and Osten W. eds Fringe 93, pp. 121-135, 1993. [90] D. Ghiglia, G. Mastin and L. Romero, ?Cellular-automata method for phase unwrapping.? Journal of the Optical Society of America A, Vol. 4, pp. 267-280. 1987. [91] A. Spik and D. Robinson, ?Investigation of the cellular automata method for phase unwrapping and its implementation on an array processor. Optics and Lasers in Engineering, Vol. 14, pp. 25-37, 1991. 166 [92] T. Flynn, ?Two-dimensional phase unwrapping with minimum weighted discontinuity,? Journal of the Optical Society of America A, Vol. 14, No. 10, pp. 2692-2701, 1997. [93] N. Ching, D. Rosenfeld and M. Braun, ?Two-dimensional phase unwrapping using a minimum spanning tree algorithm?, IEEE Trans. Image Processing, No. 1, pp. 355-365, 1992. [94] L. An, Q. Xiang, S. Chavex, ?A fast implementation of the minimum spanning tree method?, IEEE Transaction on medical imaging, Vol. 19, No. 8, pp. 805-808, 2000. [95] D. Ghiglia and L. Romero, ?Robust two-dimensinal weighted and unweighted phase unwrapping that uses fast transforms and iterative methods?, Journal of the Optical Society of America A, Vol. 11, No. 1, pp. 107-117, 1994. [96] G. Kaufmann, G. Galizzi and P. Ruiz, ?Evaluation of a preconditioned conjugate- gradient algorithm for weighted least-squares unwrapping of digital speckle- pattern interferometry phase maps,? Applied Optics, Vol. 37, No. 14, pp. 3076- 3084, 1998. [97] A. Baldi, F. Bertolino and F. Ginesu, ?On the performance of some unwrapping algorithms?, Optics and Lasers in Engineering, Vol. 37, pp. 313-330, 2002. [98] A. Anand and J. Wang, ?Strain contouring using Gabor filters: principle and algorithm?, Optical Engineering, Vol. 41, No. 6, pp. 1400-1405, 2002. [99] R. Stanford, ?Application of the least-square method to photoelastic analysis,? Experimental Mechanics, No. 2, pp. 192-197, 1980. 167 [100] A. Weeks, H. Myler and H. Wenaas, ?Computer-generated noise images for the evaluation of image processing algorithms?, Optical Engineering, Vol. 32, No. 5, pp. 982-992, 1993 [101] R. Lu, Principle and application of computer graphics, Beijing: THU press, 1996 [102] H. Taksaki, ?Moir? topography,? Applied Optics, Vol. 9, No. 6, pp. 1467-1472, 1970. [103] T. Yoshizawa and T. Tomisawa, ?Shadow moir? topography by means of the phase-shifting method,? Optical Engineering, Vol. 32, No. 7, pp. 1668-1674, 1993. [104] Y. Choi and S. Kim, ?Phase-shifting grating projection moir? topography,? Optical Engineering, Vol. 37, No. 3, pp. 1005-1010, 1998. [105] P. de Groot and L. Deck, ?Surface profiling by analysis of white-light interferograms in the spatial frequency domain,? Journal of Modern Optics, Vol. 42, No. 2, pp. 389-401, 1995. [106] L. Deck and P. de Groot, ?High-speed non-contact profiler based on scanning white light interferometry,? Applied Optics, Vol. 33, No. 31, pp. 7334-7338, 1995. [107] C. Munnerlyn and M. Latta, ?Rough surface interferometry using CO2 laser source,? Applied Optics, Vol. 7, No. 9, pp. 1858-1859, 1968. [108] O. Kwon, J. C. Wyant and C. R. Hayslett, ?Rough surface interferometry at 10.6 ?m,? Applied Optics, Vol. 19, No. 11, pp. 1862-1869, 1980. 168 [109] J. Lewandowski, B. Mongeau, M. Cormier and J. Lapierre, ?Infrared interferometers at 10 mm,? Journal of Applied Physics, Vol. 60, No. 10, pp. 3407- 3413, 1986. [110] J. Sinha and H. Tippur, ?Infrared interferometry for rough surface measurements: application to failure characterization and flaw detection,? Optical Engineering, Vol. 38, No. 8, pp. 2233-2239, 1997. [111] K. Verma and B. Han, ?Warpage measurement on dielectric rough surfaces of microelectronics devices by far infrared Fizeau interferometry,? Journal of Electronic Packaging, Transaction of the ASME, Vol. 122. No. 3, pp. 227-232, 2000. [112] K. Verma and B. Han, ?Far infrared Fizeau interferometry,? Applied Optics, Vol. 40, No. 28, pp. 4981-4987, 2001. [113] J. Blinn, ?Models of light reflection for computer synthesized pictures,? Computer Graphics, Vol. 11, No. 2, pp. 192-198, 1977. [114] R. Cook and K. Torrance ?A reflectance model for computer graphics,? ACM Transactions on Graphics, Vol. 1, No. 1, pp. 7-24, 1982. [115] X. He, K. Torrance, F. Sillion and D. Greenberg, ?A comprehensive physical model for light reflection,? Computer Graphics, Vol. 25, August 1991. [116] E. Hecht, Optics, Addison-Wesley Publishing Company, Reading, MA, 1998. [117] D. Post, B. Han and P. Ifju, High Sensitivity Moir?, Springer-Verlag, Inc., New York, 1994. 169 [118] F.-L. Dai, J. McKelvie and D. Post, "An interpretation of moir? interferometry from wavefront interference theory," Optics and Lasers in Engineering, Vol. 12, No. 2, pp. 101-118, 1990. [119] P. Hariharan, "Quasi-heterodyne Hologram Interferometry," Optical Engineering, Vol. 24, no. 4, 632-638, 1985. [120] H. Ding, R. Powell, C. Hanna and I. Ume, ?Warpage measurement comparison using shadow moir? and projection moir? methods?, IEEE Tran. On Components and Packaging Technologies, Vol. 25, No. 4, pp. 714-721, 2002. [121] B. Han and Y. Guo, ?Photomechanics tools as applied to electronic packaging product development,? in B. Han, R. Mahajan, and D. Barker (eds.), Experimental/Numerical Mechanics in Electronics Packaging, Vol. 1, pp. 11?15, Society for Experimental Mechanics, Bethel, CT, 1997. [122] M. Variyam and B. Han, ?DMD module calibration using interferometry,? TI Technical Journal, pp. 1?8, 2001. [123] Y. Guo, ?Applications of shadow moir? method in determination of thermal deformations in electronic packaging,? Proceedings of the 1995 SEM Spring Conference, Grand Rapids, MI, 1995. [124] C. Yeh, I. Ume, R. Fulton, K. Wyatt and J. Stafford, ?Correlation of analytical and experimental approaches to determine thermally induced PWB?, IEEE Trans. Components, Hybrids and Manufacturing Tech., Vol. 16, No. 8, pp. 986?995, 1993. 170 [125] B. Han, ?Thermal stresses in microelectronics subassemblies: quantitative characterization using photomechanics methods,? Journal of Thermal Stresses, Vol. 26, pp. 583-613, 2003. [126] B. Han, Y. Guo, and H.-C. Choi, ?Out-of-plane displacement measurement of printed circuit board by shadow moir? with variable sensitivity,? Proceedings of the 1993 ASME International Electronics Packaging Conference, Binghamton, NY, September 1993. [127] M. Stiteler, I. Ume and B. Leutz, ?In-process board warpage measurement in a lab scale wave soldering oven,? IEEE Trans. Components, Packaging, and Manufacturing Technology Part A, Vol. 19, No. 4, pp. 562?569, 1996. [128] H. Hertz, Zeitschrift f?r Mathematik und Physik, Vol. 28, 1883. [129] J. Mitchell, ?Elementary distributions of plane stress,? Proceedings of the London Mathematical Society, Vol. 32, pp. 35-61, 1901. [130] N. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity, P. Noordhoff, Groningen, The Netherlands, pp. 330-334, 1963. [131] I. Sokolnikoff, Mathematical Theory of Elasticity, Krieger Publishing, Malabar, pp. 283-284, 1983. [132] C. Sciammarella, Theoretical and Experimental Study on Moir? Fringes, Ph.D. Dissertation, Illinois Institute of Technology, Chicago, IL, pp. 185-192, 1960. [133] A. Durelli and L. Ferrer, ?New methods to determine elastic constants,? Materials Research and Standards, Vol. 3, pp. 988-991, 1963. 171 [134] J. C?rdenas-Garc?a, ?The moir? circular disc: two inverse problems,? Mechanics Research Communications, Vol. 28, No. 4, pp. 381-388, 2001. [135] D. Post, B. Han and P. Ifju, Chap. 7, ?Moir? methods for engineering and science ? moir? interferometry and shadow moir?,? Photomechanics for Engineers, Pramod Rastogi, ed., Springer-Verlag, 2000. [136] B. Han, D. Post and P. Ifju, ?Moir? interferometry for engineering mechanics: current practice and future development,? Journal of Strain Analysis, Vol. 36, No. 1, pp. 101-117, 2001. [137] M. Kujawinska, ?Spatial phase measurement methods,? In: Robinson D W, G T Reid, eds. Interferogram analysis. London: IOP Publishing, pp. 141~193, 1993. [138] S. Timoshenko and J. Goodier, Theory of Elasticity, Third edition, McGraw-Hill, New York, pp. 122-127, 1970. [139] R. Sanford, ?Application of the least squares method to photoelastic analysis,? Experimental Mechanics, Vol. 20, No. 6, pp. 192-197, 1980. [140] R. Sanford and J. Dally, ?A general method for determining mixed-mode stress intensity factors from isochromatic fringe patterns,? Engineering Fracture Mechanics, Vol. 11, No. 4, pp. 621-633, 1979. [141] J. Dally and W. Riley, Experimental Stress Analysis, Third edition, McGraw-Hill, New York, pp. 616-619, 1991. 172 [142] W. Gander and J. Hrebicek, Solving Problems in Scientific Computing Using Maple and MATLAB , 3rd Edition, Springer-Verlag, New York, pp. 136-139, 1997. [143] Metals Handbook, Volume 2: Properties and selection: non-ferrous alloys and special purpose materials, Tenth Edition, ASM International Materials Park, OH, 1990. [144] A. Bastawros and A. Voloshin, ?Transient thermal strain measurements in electronic packages,? IEEE Trans. Components, Hybrids and Manufacturing Technology, Vol. 13, No. 4, pp. 961-966, 1990. [145] A. Bastawros and A. Voloshin, ?In situ calibration of stress chips,? IEEE Trans. Components, Hybrids and Manufacturing Technology, Vol. 13, No. 4, pp. 888? 892, 1990. [146] Y. Guo, W. Chen and C. Lim, ?Experimental determination of thermal strains in semiconductor packaging using moir? interferometry,? Proceedings of 1992 Joint ASME/JSME Conference on Electronic Packaging, pp. 779?784, American Society of Mechanical Engineers, New York, 1992. [147] Y. Guo, W. Chen, C. Lim and C. Woychik, ?Solder ball connect (SBC) assemblies under thermal loading: I. deformation measurement via moir? interferometry, and its interpretation,? IBM Journal of Research and Development, Vol. 37, No. 5, pp. 635?648, 1993. 173 [148] T. Wu, Y. Guo and W. Chen, ?Thermal-mechanical strain characterization for printed wiring boards,? IBM Journal of Research and Development, Vol. 37, No. 5, pp. 621?634, 1993. [149] P. Tsao and A. Voloshin, ?Manufacturing stresses in the die due to die-attach process,? IEEE Trans. Components, Packaging and Manufacturing Technology Part A, Vol. 18, No. 1, pp. 201?205, 1995. [150] B. Han and Y. Guo, ?Determination of effective coefficient of thermal expansion of electronic packaging components: a whole-field approach,? IEEE Trans. Components, Packaging and Manufacturing Technology Part A, Vol. 19, No. 2, pp. 240?247, 1996. [151] B. Han, ?Deformation mechanism of two-phase solder column interconnections under highly accelerated thermal cycling condition: an experimental study,? Journal of Electronic Packaging, Trans. ASME, Vol. 119, pp. 189?196, 1997. [152] A. Voloshin, P. Tsao and R. Pearson, ?In situ evaluation of residual stresses in an organic die-attach adhesive,? Journal of Electronic Packaging, Trans. ASME, Vol. 120, No. 3, pp. 314?318, 1998. [153] S.-M. Cho, S. Cho and B. Han, ?Observing real-time thermal deformations in electronic packaging,? Experimental Techniques, Vol. 26, No. 3, pp. 25?29, 2002. [154] E. Stellrecht, B. Han and M. Pecht, ?Measurement of the hygroscopic swelling coefficient in mold compounds using moir? interferometry,? Experimental Techniques, Vol. 27, No. 4, pp. 40-44, 2003. 174 [155] S. Cho and B. Han, ?Effect of underfill on flip-chip solder bumps: an experimental study by microscopic moir? interferometry,? International Journal of Microcircuits and Electronic Packaging, Vol. 24, No. 3, pp.217?239, 2001. [156] B. Han and P. Kunthong, ?Micro-mechanical deformation analysis of surface laminar circuit in organic flip-chip package: an experimental study,? Journal of Electronic Packaging, Trans. ASME, Vol. 122, No.4, pp. 294?300, 2000. [157] B. Han, ?High sensitivity moir? interferometry for micromechanics studies,? Optical Engineering, Vol. 31, No.7, pp. 1517?1526, 1992. [158] B. Han, Z. Wu and S. Cho, ?Measurement of thermal expansion coefficient of flexible substrate,? Experimental Techniques, Vol. 25, No. 3, 22-25, 2001.