[备份]基于白光干涉实现的表面形貌仪

capric98

November 5, 2018

2658 views

No comments

5878 words

物理

基于白光干涉实现的表面形貌仪

原理

白光相干性较低、相干长度很小，使用白光的迈克尔逊干涉仪时，仅在光程差为0附近有较明显的干涉图样，呈现出双缝干涉套上高斯波包的形式
在光程差0附近的图样，可以看作是不同波数$k$在相位$\phi$下叠加得到
对于白光干涉图样的离散采样，通过傅里叶变换应可以找到一个平均的波数成分$k_0$（对所有像素按采样顺序拍成的序列分别计算FFT，然后计算k对强度的加权平均即可，注意去掉直流成分）
可以将$\phi$写成：
$$\phi=\phi_0+(k-k_0)\frac{d\phi}{d k}\bigg|_{k_0}+\frac{(k-k_0)^2}{2}\frac{d^2\phi}{d k^2}\bigg|_{k_0}+...$$
其中$\phi_0=k_0Z_0$
于是$G_0=\frac{d\phi}{dk}\bigg|_{k_0}=Z_0 + k \frac{d Z}{dk}\bigg|_{k_0}$，其中$G_0$是$Z_0$处的群速度，代回第一个式子后发现$G_0$也是$\phi$随$k$变化的一阶项系数
通过FFT，我们可以得到波数$k$及其对应的相位信息
在$k_0$附近，相位应该是随着$k$递增的，因为$\phi=kZ$，而$Z$与物体实际的高度有关（需要单独处理，因为FFT计算得到复数，使用np.angle()得到的角度在[-$\pi$,$\pi$]范围内，不保证递增性）
简单来说，通过计算$2nh=\frac{d\phi}{dk}\bigg|_{k_0}$已经可以知道物体的高度了，但有个问题是对于不同的h，$\phi$是周期性变化的，这样的算法不能区分出$h$和$h+\frac{\lambda}{2n}$，因此需要做一次Unwrap操作
论文Surface Profiling by Analysis of White-light Interferograms in the Spatial Frequency Domain给出的方法是，采用上面计算的h，考虑$\phi_0$的信息再进行一次计算（文中使用拟合的方法计算$G_0$，拟合同时给出了截距信息$\phi_0$），采用公式：
$h'=\frac{1}{2n}\{\frac{\phi_0-\alpha}{k_0}-\frac{2\pi}{k_0}Int[\frac{(\phi_0-\alpha)-(2nhk_0)}{2\pi}]\}$，其中$Int[a]$表示不超过a的最大整数
实质上是一次Unwrap操作，在python里其实有更好的二维Unwrap函数可供使用

计算部分代码

import gzip
import pickle
import os
import time
import pylab
import numpy as np

import matplotlib.pyplot as plt
# print(plt.style.available)  # show available template
plt.style.use(['seaborn-ticks','seaborn-paper'])  # use a templet

import matplotlib as mpl
# mpl.rcParams['lines.linewidth'] = 2
# mpl.rcParams['lines.color'] = 'r'
params = {
    'figure.figsize': [8, 6], # Note! figure unit is inch!  scale fontz size 2.54 to looks like unit cm
    'axes.labelsize': 7*2.54, # scale 2.54 to change to figure unit looks as cm
    'font.size':  7*2.54,
    'lines.linewidth':2,
    'legend.fontsize': 7*2.54,
    'xtick.labelsize': 6*2.54,
    'ytick.labelsize': 6*2.54,
    'text.usetex': False,  
    'xtick.direction': "in",
    'ytick.direction': "in", # ticket inside
    'legend.frameon' : True, 
    'legend.edgecolor': 'black',
    'legend.shadow': True,
    'legend.framealpha':1,
#     'patch.linewidth' : 0.5, 
}
mpl.rcParams.update(params)

with gzip.open('data.pklz', 'rb') as f:
    frames = pickle.load(f)

frames_N = len(frames)

Pi = np.pi

def phase_unwrap(Raw):
    Len_Raw = len(Raw)
    Threshold = 1.5
    Unwraped = Raw
    for i in range(Len_Raw-1):
        Unwraped[i+1] = Raw[i+1] + 2*Pi*np.rint((Unwraped[i]-Raw[i+1])/(2.0*Pi))
        if (Unwraped[i+1]<Unwraped[i])and(np.abs(Unwraped[i]-Unwraped[i+1])>Threshold) : 
          Unwraped[i+1] += 2*Pi
    return Unwraped

def Smooth(Raw):
    Threshold = 0.5
    Smooth_Img = Raw
    Width , Height = Raw.shape
    for i in range(Width-2):
        for j in range(Height-2):
            Tmp = (Raw[i,j+1] + Raw[i+2,j+1] + Raw[i+1,j] + Raw[i+1,j+2])/4.0
            if np.abs(Tmp-Smooth_Img[i+1,j+1])>Threshold :
                Smooth_Img[i+1,j+1] = Tmp
    return Smooth_Img

def ProfileCalculation(frames, Step, fit_window=2, SkipZero=5, Threshold = 1/10.0):

    Width , Height = frames[0].shape
    
    N = len(frames)
    Len = Width * Height

    lines = [frame.reshape(-1, 1) for frame in frames]
    lines = np.concatenate(tuple(lines), axis=1)
    
    Angles = []
    Data = []  #Store the phase & k after FT
    
    Bar = 0
    Count =0
    Part = 0 #Relate to Progress Bar
    
    for i, line in enumerate(lines):
        line = np.concatenate((np.ones(len(line),dtype=np.uint8)*int(np.average(line)),line))
        rfft = np.fft.rfft(line)
        rfft_abs = np.abs(rfft)
        kn = np.argmax(rfft_abs[SkipZero:-1]) + SkipZero 
        Data.append([rfft_abs[kn], kn])
        Angles.append(np.angle(rfft)) # angle is in range(-Pi,Pi)
        Count += 1
        if int(100.0*Count / Len)!=Bar :
            Bar +=1
            #os.system("cls")
            print("Processing ... %.2f%%" %(50*(Count / Len) + Part)) 
    Part = 50
    Bar =0
    Count =0

    Data = np.array(Data)

    kn_Max_Value = np.max(Data[:, 0])
    Mask = Data[:, 0] > kn_Max_Value * Threshold  
        #Throw the background data which the value(k=kn) less than max_value /10.0
    Data = Data[Mask, :]

    k0 = np.sum(Data[:, 0]*Data[:, 1])/np.sum(Data[:, 0]) #Calculate the most efficient k0 of white light
    k0_int = int(k0) 
    Fit_Range = np.array(range((k0_int - fit_window), (k0_int + fit_window + 1)))

    Image = []

    for i, angle in enumerate(Angles):
        if Mask[i] :
            Count += 1
            angle = phase_unwrap(angle[Fit_Range])
                #Using np.angle we get the phase in range of (-Pi,Pi)
                #Around the k0, the phi is equal k times Z,
                #since the Z is a const relate to the real distance of the object
                #the phase should be increasing.
            k , Slope , Phi0 = (np.polyfit(Fit_Range-k0,angle,2)) 
                #Polyfit around k0
                #Slope -> G0
                #Fit_Range-k0 so that intercept is Phi0 in k=k0
            Image.append(Phi0-2*np.pi*np.rint((Phi0-k0_int * Slope)/(2*np.pi)))
        else :
            Count += 1
            Image.append(0)
                #Backgroud is simply set as h=0
        if int(100*Count /Len)!=Bar :
            Bar += 1
            #os.system("cls")
            print("Processing ... %.2f%%" %(50*(Count / Len)+ Part))

    Image = np.array(Image).reshape(Width,Height)
    #for i in range(Width):
    #    Image[i,:] = unwrap_phase(Image[i,:])
    Image = unwrap_phase(Image)
    k = 2 * np.pi * k0 / (Step * N)
    Image = Image/k/2
    return Image

from skimage.restoration import unwrap_phase
    
image = ProfileCalculation(frames[:], Step=0.1)
image = Smooth(image)

img = plt.imshow(image)
img.set_cmap('Spectral_r')
plt.colorbar()
pylab.show()

[备份]基于白光干涉实现的表面形貌仪