優化星形輪廓（高斯）擬合算法

Question

我需要將幾千個二維高斯函數擬合到CCD圖像上的星形輪廓（14x14像素塊）上，並獲得質心座標，長軸和短軸的FWHM以及長軸的旋轉角度。問題是我的當前代碼執行時間太長。在i7處理器上有幾十秒的時間，我需要它做得更快。儘可能快。我測試了幾個高斯擬合函數，看起來在AsPyLib中使用的函數是最快的http://www.aspylib.com/doc/aspylib_fitting.html 下面是我試圖使運行速度更快的代碼。分析顯示大部分時間都花費在mplfit函數中。所以我的問題是如果這可以加速嗎？我嘗試了對代碼進行cython化，但它提供了非常小的提升。計算時刻（快10倍）不適用於我有許多圖像，由於噪聲問題，使估計不可靠。可能由於星形輪廓由於像差而經常遠離高斯。 多處理不是一個解決方案，或者是因爲新的流程創建開銷太高而不適合僅僅一個星形配置文件。 任何想法在哪裏看得更遠？

import numpy as np 
from scipy.optimize import leastsq 

def fit_gauss_elliptical(xy, data): 
    """ 
    --------------------- 
    Purpose 
    Fitting a star with a 2D elliptical gaussian PSF. 
    --------------------- 
    Inputs 
    * xy (list) = list with the form [x,y] where x and y are the integer positions in the complete image of the first pixel (the one with x=0 and y=0) of the small subimage that is used for fitting. 
    * data (2D Numpy array) = small subimage, obtained from the full FITS image by slicing. It must contain a single object : the star to be fitted, placed approximately at the center. 
    --------------------- 
    Output (list) = list with 8 elements, in the form [maxi, floor, height, mean_x, mean_y, fwhm_small, fwhm_large, angle]. The list elements are respectively: 
    - maxi is the value of the star maximum signal, 
    - floor is the level of the sky background (fit result), 
    - height is the PSF amplitude (fit result), 
    - mean_x and mean_y are the star centroid x and y positions, on the full image (fit results), 
    - fwhm_small is the smallest full width half maximum of the elliptical gaussian PSF (fit result) in pixels 
    - fwhm_large is the largest full width half maximum of the elliptical gaussian PSF (fit result) in pixels 
    - angle is the angular direction of the largest fwhm, measured clockwise starting from the vertical direction (fit result) and expressed in degrees. The direction of the smallest fwhm is obtained by adding 90 deg to angle. 
    --------------------- 
    """ 

    #find starting values 
    dat=data.flatten() 
    maxi = data.max() 
    floor = np.ma.median(dat) 
    height = maxi - floor 
    if height==0.0:    #if star is saturated it could be that median value is 32767 or 65535 --> height=0 
     floor = np.mean(dat) 
     height = maxi - floor 

    mean_x = (np.shape(data)[0]-1)/2 
    mean_y = (np.shape(data)[1]-1)/2 

    fwhm = np.sqrt(np.sum((data>floor+height/2.).flatten())) 
    fwhm_1 = fwhm 
    fwhm_2 = fwhm 
    sig_1 = fwhm_1/(2.*np.sqrt(2.*np.log(2.))) 
    sig_2 = fwhm_2/(2.*np.sqrt(2.*np.log(2.)))  

    angle = 0. 

    p0 = floor, height, mean_x, mean_y, sig_1, sig_2, angle 

    #--------------------------------------------------------------------------------- 
    #fitting gaussian 
    def gauss(floor, height, mean_x, mean_y, sig_1, sig_2, angle): 

     A = (np.cos(angle)/sig_1)**2. + (np.sin(angle)/sig_2)**2. 
     B = (np.sin(angle)/sig_1)**2. + (np.cos(angle)/sig_2)**2. 
     C = 2.0*np.sin(angle)*np.cos(angle)*(1./(sig_1**2.)-1./(sig_2**2.)) 

     #do not forget factor 0.5 in exp(-0.5*r**2./sig**2.)  
     return lambda x,y: floor + height*np.exp(-0.5*(A*((x-mean_x)**2)+B*((y-mean_y)**2)+C*(x-mean_x)*(y-mean_y))) 

    def err(p,data): 
     return np.ravel(gauss(*p)(*np.indices(data.shape))-data) 

    p = leastsq(err, p0, args=(data), maxfev=200) 
    p = p[0] 

    #--------------------------------------------------------------------------------- 
    #formatting results 
    floor = p[0] 
    height = p[1] 
    mean_x = p[2] + xy[0] 
    mean_y = p[3] + xy[1] 

    #angle gives the direction of the p[4]=sig_1 axis, starting from x (vertical) axis, clockwise in direction of y (horizontal) axis 
    if np.abs(p[4])>np.abs(p[5]): 

     fwhm_large = np.abs(p[4]) * (2.*np.sqrt(2.*np.log(2.))) 
     fwhm_small = np.abs(p[5]) * (2.*np.sqrt(2.*np.log(2.))) 
     angle = np.arctan(np.tan(p[6])) 

    else: #then sig_1 is the smallest : we want angle to point to sig_y, the largest 

     fwhm_large = np.abs(p[5]) * (2.*np.sqrt(2.*np.log(2.))) 
     fwhm_small = np.abs(p[4]) * (2.*np.sqrt(2.*np.log(2.))) 
     angle = np.arctan(np.tan(p[6]+np.pi/2.)) 

    output = [maxi, floor, height, mean_x, mean_y, fwhm_small, fwhm_large, angle] 
    return output 

Answer 1

我也碰到過這種方法之前，但沒有實現它：

http://www.nature.com/nmeth/journal/v9/n7/full/nmeth.2071.html

可能是值得一試的快捷方法。甚至有代碼複製到他的網站上http://physics-server.uoregon.edu/~raghu/particle_tracking.html