4

BACKGROUND I am writing a program to calculate surface brightness of a galaxy as a function of elliptical radius. This first involves reading in a .fits file, which is stored in a numpy array such that array[x][y] will return the surface brightness value at that (x,y) pixel.

To calculate surface brightness, I need to be able to fit an ellipse to the galaxy at some minimum size and find the median surface brightness within that ellipse, then increase the size of the ellipse and find the surface brightness for each annuli. Different sized annuli will be looped through until the surface brightness drops beneath a certain ratio to the background noise.

PROBLEM Given parameters for the ellipse include position angle, center coordinates for the x and y pixels, and ratio of B/A. I am having trouble finding any sort of method which will allow me to fit an ellipse to essentially an array of arrays. Help please??

Garrett Fitzgerald
  • 49
  • 3
  • What kind of trouble do you have? – moffeltje Jul 01 '15 at 13:46
  • [This](http://scikit-image.org/docs/dev/auto_examples/plot_circular_elliptical_hough_transform.html) may help. – Jaime Jul 01 '15 at 15:57
  • It sounds like you're doing multi-gaussian expansion (http://www-astro.physics.ox.ac.uk/~mxc/software/) or something similar. If this is not what you are looking for, please clarify "fit an ellipse to an array of arrays". – keflavich Jul 02 '15 at 12:14
  • Okay so "fit an ellipse to an array of arrays" is a poor explanation. I am writing a program that will complete aperture photometry on a galaxy given the center coordinates, the position angle of the ellipse, and the ratio of the ellipse in terms of B/A. Basically, I need to be able to include the pixels within the ellipse and find the median value out of those pixels- this is the surface brightness at that given radius. Then increase the radius of the ellipse and find the surface brightness at each radius. Output is flux/area (surface brightness) at each elliptical annuli (radius). Thank you. – Garrett Fitzgerald Jul 02 '15 at 15:11

1 Answers1

1

I think I may have worked on a problem similar to yours in a research project of mine involving Laue Diffraction patterns (material science). I was tasked with finding the length, width, and angle of inclination of each peak in a diffraction pattern given the center coordinates of each peak. My solution was to select a region of interest around the peak and threshold, filter, etc such that there was only one peak in the sub-image. I then built a function to fit these parameters to the resulting ellipse:

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

I just threw this together in case it was similar to what you were looking for. If you need more explanation, feel free to comment. The sources I used (math) are within the comments.

Aaron
  • 10,133
  • 1
  • 24
  • 40