4

背景 我正在编写一个程序来计算星系的表面亮度作为椭圆半径的函数。这首先涉及读取 .fits 文件，该文件存储在一个 numpy 数组中，这样 array[x][y] 将返回该 (x,y) 像素的表面亮度值。

为了计算表面亮度，我需要能够以某个最小尺寸将椭圆拟合到星系中，并找到该椭圆内的表面亮度中值，然后增加椭圆的大小并找到每个环的表面亮度。不同大小的环将循环通过，直到表面亮度下降到与背景噪声的某个比率以下。

问题 椭圆的给定参数包括位置角、x 和 y 像素的中心坐标以及 B/A 比。我很难找到任何可以让我将椭圆拟合到本质上的数组数组的方法。请帮忙？？

4

1 回答 1

1

我想我可能在我的一个涉及劳厄衍射模式（材料科学）的研究项目中解决了与您类似的问题。在给定每个峰的中心坐标的情况下，我的任务是在衍射图中找到每个峰的长度、宽度和倾斜角度。我的解决方案是在峰值和阈值、过滤器等周围选择一个感兴趣的区域，这样子图像中只有一个峰值。然后我构建了一个函数来将这些参数拟合到生成的椭圆中：

from xml.etree.cElementTree import parse
import numpy as np
from os import listdir, getcwd
from scipy import ndimage
import h5py
import multiprocessing
import time
#from dicttoxml import dicttoxml
#from xml.dom.minidom import parseString
import sys
import cPickle as pickle
from threading import Thread
from skimage.measure import moments

def fitEllipse(data):
    '''
    Returns the length of the long and short axis and the angle measure
    of the long axis to the horizontal of the best fit ellipsebased on
    image moments.

    usage: longAxis, shortAxis, angle = fitEllipse(N_by_M_image_as_array)
    '''
    # source:
    #     Kieran F. Mulchrone, Kingshuk Roy Choudhury,
    # Fitting an ellipse to an arbitrary shape:
    # implications for strain analysis, Journal of
    # Structural Geology, Volume 26, Issue 1,
    # January 2004, Pages 143-153, ISSN 0191-8141,
    # <http://dx.doi.org/10.1016/S0191-8141(03)00093-2.>
    #     Lourena Rocha, Luiz Velho, Paulo Cezar P. Carvalho
    # Image Moments-Based Structuring and Tracking of
    # Objects, IMPA-Instituto Nacional de Matematica Pura
    # e Aplicada. Estrada Dona Castorina, 110, 22460
    # Rio de Janeiro, RJ, Brasil,
    # <http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1167130>

    m = moments(data, 2) # super fast compated to anything in pure python
    xc = m[1,0] / m[0,0]
    yc = m[0,1] / m[0,0]
    a = (m[2,0] / m[0,0]) - (xc**2)
    b = 2 * ((m[1,1] / m[0,0]) - (xc * yc))
    c = (m[0,2] / m[0,0]) - (yc**2)
    theta = .5 * (np.arctan2(b, (a - c)))
    w = np.sqrt(6 * (a + c - np.sqrt(b**2 + (a-c)**2)))
    l = np.sqrt(6 * (a + c + np.sqrt(b**2 + (a-c)**2)))
    return l, w, theta

我只是把它放在一起，以防它与你正在寻找的相似。如果您需要更多解释，请随时发表评论。我使用的来源（数学）在评论中。

于 2015-07-01T14:05:50.830 回答