Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fast build ellipse model #1288

Open
wants to merge 6 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
118 changes: 44 additions & 74 deletions photutils/isophote/model.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,16 +6,15 @@


import numpy as np

from .geometry import EllipseGeometry
import ctypes as ct
import os

__all__ = ['build_ellipse_model']


def build_ellipse_model(shape, isolist, fill=0., high_harmonics=False):
def build_ellipse_model(shape, isolist, nthreads=1, fill=0., high_harmonics=False):
"""
Build a model elliptical galaxy image from a list of isophotes.

For each ellipse in the input isophote list the algorithm fills the
output image array with the corresponding isophotal intensity.
Pixels in the output array are in general only partially covered by
Expand All @@ -24,27 +23,24 @@ def build_ellipse_model(shape, isolist, fill=0., high_harmonics=False):
each pixel by storing the partial area information in an auxiliary
array. The information in this array is then used to normalize the
pixel intensities.

Parameters
----------
shape : 2-tuple
The (ny, nx) shape of the array used to generate the input
``isolist``.

isolist : `~photutils.isophote.IsophoteList` instance
The isophote list created by the `~photutils.isophote.Ellipse`
class.

nthreads: float, optiomal

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
nthreads: float, optiomal
nthreads: float, optional

Number of threads to perform work. Default is 1 (serial code).

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Number of threads to perform work. Default is 1 (serial code).
Number of threads to execute the task. Default is 1 (serial code).

fill : float, optional
The constant value to fill empty pixels. If an output pixel has
no contribution from any isophote, it will be assigned this
value. The default is 0.

high_harmonics : bool, optional
Whether to add the higher-order harmonics (i.e., ``a3``, ``b3``,
``a4``, and ``b4``; see `~photutils.isophote.Isophote` for
details) to the result.

Returns
-------
result : 2D `~numpy.ndarray`
Expand Down Expand Up @@ -91,71 +87,42 @@ def build_ellipse_model(shape, isolist, fill=0., high_harmonics=False):

# correct deviations cased by fluctuations in spline solution
eps_array[np.where(eps_array < 0.)] = 0.

result = np.zeros(shape=shape)
weight = np.zeros(shape=shape)

eps_array[np.where(eps_array < 0.)] = 0.05

# for each interpolated isophote, generate intensity values on the
# output image array
# for index in range(len(finely_spaced_sma)):
for index in range(1, len(finely_spaced_sma)):
sma0 = finely_spaced_sma[index]
eps = eps_array[index]
pa = pa_array[index]
x0 = x0_array[index]
y0 = y0_array[index]
geometry = EllipseGeometry(x0, y0, sma0, eps, pa)

intens = intens_array[index]

# scan angles. Need to go a bit beyond full circle to ensure
# full coverage.
r = sma0
phi = 0.
while phi <= 2*np.pi + geometry._phi_min:
# we might want to add the third and fourth harmonics
# to the basic isophotal intensity.
harm = 0.
if high_harmonics:
harm = (a3_array[index] * np.sin(3.*phi) +
b3_array[index] * np.cos(3.*phi) +
a4_array[index] * np.sin(4.*phi) +
b4_array[index] * np.cos(4.*phi)) / 4.

# get image coordinates of (r, phi) pixel
x = r * np.cos(phi + pa) + x0
y = r * np.sin(phi + pa) + y0
i = int(x)
j = int(y)

if (i > 0 and i < shape[1] - 1 and j > 0 and j < shape[0] - 1):
# get fractional deviations relative to target array
fx = x - float(i)
fy = y - float(j)

# add up the isophote contribution to the overlapping pixels
result[j, i] += (intens + harm) * (1. - fy) * (1. - fx)
result[j, i + 1] += (intens + harm) * (1. - fy) * fx
result[j + 1, i] += (intens + harm) * fy * (1. - fx)
result[j + 1, i + 1] += (intens + harm) * fy * fx

# add up the fractional area contribution to the
# overlapping pixels
weight[j, i] += (1. - fy) * (1. - fx)
weight[j, i + 1] += (1. - fy) * fx
weight[j + 1, i] += fy * (1. - fx)
weight[j + 1, i + 1] += fy * fx

# step towards next pixel on ellipse
phi = max((phi + 0.75 / r), geometry._phi_min)
r = max(geometry.radius(phi), 0.5)
# if outside image boundaries, ignore.
else:
break

# zero weight values must be set to 1.

# convert everything to C-type array (pointers)
c_fss_array = ct.c_void_p(finely_spaced_sma.ctypes.data)
c_intens_array = ct.c_void_p(intens_array.ctypes.data)
c_eps_array = ct.c_void_p(eps_array.ctypes.data)
c_pa_array = ct.c_void_p(pa_array.ctypes.data)
c_x0_array = ct.c_void_p(x0_array.ctypes.data)
c_y0_array = ct.c_void_p(y0_array.ctypes.data)
c_a3_array = ct.c_void_p(a3_array.ctypes.data)
c_b3_array = ct.c_void_p(b3_array.ctypes.data)
c_a4_array = ct.c_void_p(a4_array.ctypes.data)
c_b4_array = ct.c_void_p(b4_array.ctypes.data)

# initialize result and weight arrays, also as 1D ctype array
result = np.zeros(shape=(shape[1]*shape[0],))
weight = np.zeros(shape=(shape[1]*shape[0],))
c_result = ct.c_void_p(result.ctypes.data)
c_weight = ct.c_void_p(weight.ctypes.data)

# convert high_harmnics bool flag to int,
# convert all other ints to ctype
c_high_harm = ct.c_int(int(high_harmonics))
c_N = ct.c_int(len(finely_spaced_sma))
c_nrows = ct.c_int(shape[0])
c_ncols = ct.c_int(shape[1])

# load into C worker function (worker.so should be in same directory)
lib = ct.cdll.LoadLibrary(os.path.dirname(os.path.abspath(__file__)) + '/worker.so')
lib.worker.restype = None
lib.worker(c_result, c_weight, c_nrows, c_ncols, c_N, c_high_harm,
c_fss_array, c_intens_array, c_eps_array, c_pa_array,
c_x0_array, c_y0_array, c_a3_array, c_b3_array,
c_a4_array, c_b4_array, ct.c_int(nthreads))

# zero weight values must be set to 1.
weight[np.where(weight <= 0.)] = 1.

# normalize
Expand All @@ -164,4 +131,7 @@ def build_ellipse_model(shape, isolist, fill=0., high_harmonics=False):
# fill value
result[np.where(result == 0.)] = fill

return result
# reshape
result = result.reshape(shape[0], shape[1])

return result
47 changes: 47 additions & 0 deletions photutils/isophote/tests/test_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
from astropy.utils.data import get_pkg_data_filename
import numpy as np
import pytest
import multiprocessing as mp

from .make_test_data import make_test_image
from ..ellipse import Ellipse
Expand Down Expand Up @@ -43,6 +44,34 @@ def test_model():
residual = data - model
assert np.mean(residual) <= 5.0
assert np.mean(residual) >= -5.0

@pytest.mark.remote_data
@pytest.mark.skipif('not HAS_SCIPY')
def test_model_parallel():
path = get_path('isophote/M105-S001-RGB.fits',
location='photutils-datasets', cache=True)
hdu = fits.open(path)
data = hdu[0].data[0]
hdu.close()

g = EllipseGeometry(530., 511, 10., 0.1, 10./180.*np.pi)
ellipse = Ellipse(data, geometry=g, threshold=1.e5)
# NOTE: this sometimes emits warnings (e.g., py38, ubuntu), but
# sometimes not. Here we simply ignore any RuntimeWarning, whether
# there is one or not.
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
isophote_list = ellipse.fit_image()

# test parallel on a thread pool, w/ size equal to CPU thread count
model = build_ellipse_model(data.shape, isophote_list, nthreads=mp.cpu_count(),
fill=np.mean(data[10:100, 10:100]))

assert data.shape == model.shape

residual = data - model
assert np.mean(residual) <= 5.0
assert np.mean(residual) >= -5.0


@pytest.mark.skipif('not HAS_SCIPY')
Expand All @@ -62,6 +91,24 @@ def test_model_simulated_data():
assert np.mean(residual) <= 5.0
assert np.mean(residual) >= -5.0

@pytest.mark.skipif('not HAS_SCIPY')
def test_model_large_array():
data = make_test_image(nx=5000, ny=5000, i0=100., sma=500., eps=0.3,
pa=np.pi/4., noise=0.05, seed=0)

g = EllipseGeometry(2500., 2500., 100., 0.3, np.pi/4.)
ellipse = Ellipse(data, geometry=g, threshold=1.e5)
isophote_list = ellipse.fit_image()
# test parallel on a thread pool, w/ size equal to CPU thread count
model = build_ellipse_model(data.shape, isophote_list, nthreads=mp.cpu_count(),
fill=np.mean(data[0:50, 0:50]))

assert data.shape == model.shape

residual = data - model

assert np.mean(residual) <= 5.0
assert np.mean(residual) >= -5.0

@pytest.mark.skipif('not HAS_SCIPY')
def test_model_minimum_radius():
Expand Down
Binary file added photutils/isophote/worker.so
Binary file not shown.