Skip to content

Bug fixes #166

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.

Sign up for GitHub

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

Merged
PJ-Watson merged 2 commits into from
Feb 18, 2025
Merged

Bug fixes #166

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0b1f4e5
Update benchmark script, expose arguments.
PJ-Watson Feb 18, 2025
dfa868b
Add descriptive warning for aperture photometry.
PJ-Watson Feb 18, 2025
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
152 changes: 102 additions & 50 deletions bench.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
from __future__ import print_function

import time
from os.path import join
from pathlib import Path

import numpy as np

Expand All @@ -30,44 +30,99 @@
from astropy.io.fits import getdata

HAVE_FITS = True
NEED_BYTESWAP = True
except ImportError:
HAVE_FITS = False

CONDENSED = True
import argparse

parser = argparse.ArgumentParser(description="SEP Benchmark script.")
parser.add_argument(
"-c",
"--condensed",
default=False,
action="store_true",
help="Only run a condensed subset of the benchmarks.",
)
parser.add_argument(
"-t",
"--tiles",
default=4,
type=int,
help=(
"The maximum number of tiles for the image background benchmark, "
"i.e. the image size will be `(t*256)^2`."
),
)
parser.add_argument(
"-a",
"--apertures",
default=1000,
type=int,
help=(
"The number of apertures to use for benchmarking the aperture "
"photometry, by default 1000."
),
)
parser.add_argument(
"-r",
"--radii",
default=[3, 5, 10, 20],
type=float,
nargs="+",
help=(
"The radii of the apertures to test. Ignored if the `--condensed` "
"flag is passed, otherwise defaults to `3, 5, 10, 20`."
),
)
parser.add_argument(
"-n",
"--nloop",
default=50,
type=int,
help=("The number of loops to run for all tests, by default 50."),
)
args = parser.parse_args()

if HAVE_FITS:
rawdata = getdata(join("data", "image.fits")) # original is 256 x 256
data = np.tile(rawdata, (4, 4))
rawdata = getdata(
Path(__file__).parent / "data" / "image.fits"
) # original is 256 x 256
if NEED_BYTESWAP:
rawdata = rawdata.astype(rawdata.dtype.newbyteorder("="))
else:
print("No FITS reader found, generating random data instead.")
rawdata = np.random.random((256, 256))

print("test image shape:", data.shape)
print("test image dtype:", data.dtype)
data = np.tile(rawdata, (4, 4))

t0 = time.time()
bkg = sep.Background(data) # estimate background
t1 = time.time()
print("measure background: {0:6.2f} ms".format((t1 - t0) * 1.0e3))
print("test image shape:", data.shape)
print("test image dtype:", data.dtype)

t0 = time.time()
bkg.subfrom(data) # subtract it
t1 = time.time()
print("subtract background: {0:6.2f} ms".format((t1 - t0) * 1.0e3))
t0 = time.time()
bkg = sep.Background(data) # estimate background
t1 = time.time()
print("measure background: {0:6.2f} ms".format((t1 - t0) * 1.0e3))

t0 = time.time()
backarr = bkg.back(dtype=np.float64) # background
t1 = time.time()
print("background array: {0:6.2f} ms".format((t1 - t0) * 1.0e3))
t0 = time.time()
bkg.subfrom(data) # subtract it
t1 = time.time()
print("subtract background: {0:6.2f} ms".format((t1 - t0) * 1.0e3))

t0 = time.time()
rmsarr = bkg.rms()
t1 = time.time()
print("rms array: {0:6.2f} ms".format((t1 - t0) * 1.0e3))
t0 = time.time()
backarr = bkg.back(dtype=np.float64) # background
t1 = time.time()
print("background array: {0:6.2f} ms".format((t1 - t0) * 1.0e3))

t0 = time.time()
objects = sep.extract(data, 1.5 * bkg.globalrms)
t1 = time.time()
print(
"extract: {0:6.2f} ms [{1:d} objects]".format((t1 - t0) * 1.0e3, len(objects))
)
t0 = time.time()
rmsarr = bkg.rms()
t1 = time.time()
print("rms array: {0:6.2f} ms".format((t1 - t0) * 1.0e3))

t0 = time.time()
objects = sep.extract(data, 1.5 * bkg.globalrms)
t1 = time.time()
print("extract: {0:6.2f} ms [{1:d} objects]".format((t1 - t0) * 1.0e3, len(objects)))

# --------------------------------------------------------------------------
# Background subtraction
Expand All @@ -94,51 +149,48 @@
)
blankline = "| | |"

nloop = 50

for ntile in [4]:
for ntile in np.arange(1, args.tiles + 1, dtype=int):
data = np.tile(rawdata, (ntile, ntile))
line = "| {0:4d}^2 image background |".format(data.shape[0])

t0 = time.time()
for _ in range(0, nloop):
for _ in range(0, args.nloop):
bkg = sep.Background(data)
t1 = time.time()
t_sep = (t1 - t0) * 1.0e3 / nloop
t_sep = (t1 - t0) * 1.0e3 / args.nloop
line += " {0:7.2f} ms |".format(t_sep)

if HAVE_PHOTUTILS:
t0 = time.time()
for _ in range(0, nloop):
for _ in range(0, args.nloop):
from photutils import background

bkg = background.Background2D(data, (64, 64)) # estimate background
t1 = time.time()
t_pu = (t1 - t0) * 1.0e3 / nloop
t_pu = (t1 - t0) * 1.0e3 / args.nloop
line += " {0:7.2f} ms | {1:6.2f} |".format(t_pu, t_pu / t_sep)

print(line)

# ------------------------------------------------------------------------------
# Circular aperture photometry benchmarks

if not CONDENSED:
if not args.condensed:
print(blankline)
line = "| **aperture photometry** | |"
if HAVE_PHOTUTILS:
line += " | |"
print(line)

naper = 1000
data = np.ones((2000, 2000), dtype=np.float32)
x = np.random.uniform(200.0, 1800.0, naper)
y = np.random.uniform(200.0, 1800.0, naper)
x = np.random.uniform(200.0, 1800.0, args.apertures)
y = np.random.uniform(200.0, 1800.0, args.apertures)

if CONDENSED:
if args.condensed:
r_list = [5.0]
subpix_list = [(5, "subpixel", "subpix=5"), (0, "exact", "exact")]
else:
r_list = [3.0, 5.0, 10.0, 20.0]
r_list = args.radii
subpix_list = [
(1, "center", "subpix=1"),
(5, "subpixel", "subpix=5"),
Expand All @@ -151,28 +203,28 @@
line = "| circles r={0:2d} {1:8s} |".format(int(r), label)

t0 = time.time()
for _ in range(0, nloop):
for _ in range(0, args.nloop):
flux, fluxerr, flag = sep.sum_circle(data, x, y, r, subpix=subpix)
t1 = time.time()
t_sep = (t1 - t0) * 1.0e6 / naper / nloop
t_sep = (t1 - t0) * 1.0e6 / args.apertures / args.nloop
line += " {0:7.2f} us/aper |".format(t_sep)

if HAVE_PHOTUTILS:
from photutils import aperture

apertures = photutils.aperture.CircularAperture(np.column_stack((x, y)), r)
t0 = time.time()
for _ in range(0, nloop):
for _ in range(0, args.nloop):
res = photutils.aperture.aperture_photometry(
data, apertures, method=method, subpixels=subpix
)
t1 = time.time()
t_pu = (t1 - t0) * 1.0e6 / naper / nloop
t_pu = (t1 - t0) * 1.0e6 / args.apertures / args.nloop
line += " {0:7.2f} us/aper | {1:6.2f} |".format(t_pu, t_pu / t_sep)

print(line)

if not CONDENSED:
if not args.condensed:
print(blankline)

a = 1.0
Expand All @@ -184,12 +236,12 @@
line = "| ellipses r={0:2d} {1:8s} |".format(int(r), label)

t0 = time.time()
for _ in range(0, nloop):
for _ in range(0, args.nloop):
flux, fluxerr, flag = sep.sum_ellipse(
data, x, y, a, b, theta, r, subpix=subpix
)
t1 = time.time()
t_sep = (t1 - t0) * 1.0e6 / naper / nloop
t_sep = (t1 - t0) * 1.0e6 / args.apertures / args.nloop
line += " {0:7.2f} us/aper |".format(t_sep)

if HAVE_PHOTUTILS:
Expand All @@ -199,12 +251,12 @@
np.column_stack((x, y)), a * r, b * r, theta
)
t0 = time.time()
for _ in range(0, nloop):
for _ in range(0, args.nloop):
res = photutils.aperture.aperture_photometry(
data, apertures, method=method, subpixels=subpix
)
t1 = time.time()
t_pu = (t1 - t0) * 1.0e6 / naper / nloop
t_pu = (t1 - t0) * 1.0e6 / args.apertures / args.nloop
line += " {0:7.2f} us/aper | {1:6.2f} |".format(t_pu, t_pu / t_sep)

print(line)
7 changes: 7 additions & 0 deletions docs/apertures.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,12 @@ aperture photometry if the Kron radius is too small. For example:
fluxerr[use_circle] = cfluxerr
flag[use_circle] = cflag

.. warning::
Caution should be used when calculating Kron radii in crowded fields.
In almost all cases, one would want to pass in a segmentation map to
mask out nearby objects, as described below in
:ref:`segmentation masking`.

Equivalent of FLUX_RADIUS in Source Extractor
---------------------------------------------

Expand Down Expand Up @@ -179,6 +185,7 @@ of 0.5 and a normalizing flux of ``FLUX_AUTO``. The equivalent here is:
sig = 2. / 2.35 * r # r from sep.flux_radius() above, with fluxfrac = 0.5
xwin, ywin, flag = sep.winpos(data, objs['x'], objs['y'], sig)

.. _segmentation masking:

Segmentation-masked image measurements
--------------------------------------
Expand Down
30 changes: 21 additions & 9 deletions sep.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -221,10 +221,10 @@ cdef int _get_sep_dtype(dtype) except -1:
"""Convert a numpy dtype to the corresponding SEP dtype integer code."""
if not dtype.isnative:
raise ValueError(
"Input array with dtype '{0}' has non-native byte order. "
"Input array with dtype `{0}` has non-native byte order. "
"Only native byte order arrays are supported. "
"To change the byte order of the array 'data', do "
"'data = data.byteswap().newbyteorder()'".format(dtype))
"To change the byte order of the array `data`, do "
"`data = data.astype(data.dtype.newbyteorder('='))`".format(dtype))
t = dtype.type
if t is np.single:
return SEP_TFLOAT
Expand Down Expand Up @@ -1060,10 +1060,17 @@ def sum_circle(np.ndarray data not None, x, y, r,
1, SEP_MASK_IGNORE, &bkgflux, &bkgfluxerr, &bkgarea, &bkgflag)
_assert_ok(status)

if not bkgarea > 0:
raise ValueError(
"The background annulus does not contain any valid pixels, "
"for the object at index "
f"{np.PyArray_MultiIter_INDEX(it)}."
)

if area1 > 0:
flux1 -= bkgflux / bkgarea * area1
bkgfluxerr = bkgfluxerr / bkgarea * area1
fluxerr1 = sqrt(fluxerr1*fluxerr1 + bkgfluxerr*bkgfluxerr)
flux1 -= bkgflux / bkgarea * area1
bkgfluxerr = bkgfluxerr / bkgarea * area1
fluxerr1 = sqrt(fluxerr1*fluxerr1 + bkgfluxerr*bkgfluxerr)
(<double*>np.PyArray_MultiIter_DATA(it, 6))[0] = flux1
(<double*>np.PyArray_MultiIter_DATA(it, 7))[0] = fluxerr1
(<short*>np.PyArray_MultiIter_DATA(it, 8))[0] = flag1
Expand Down Expand Up @@ -1340,7 +1347,6 @@ def sum_ellipse(np.ndarray data not None, x, y, a, b, theta, r=1.0,
sumerr = np.empty(shape, dt)
flag = np.empty(shape, np.short)

#it = np.broadcast(x, y, a, b, theta, r, sum, sumerr, flag)
it = np.broadcast(x, y, a, b, theta, r, seg_id, sum, sumerr, flag)
while np.PyArray_MultiIter_NOTDONE(it):
status = sep_sum_ellipse(
Expand Down Expand Up @@ -1376,7 +1382,6 @@ def sum_ellipse(np.ndarray data not None, x, y, a, b, theta, r=1.0,
sumerr = np.empty(shape, dt)
flag = np.empty(shape, np.short)

# it = np.broadcast(x, y, a, b, theta, r, rin, rout, sum, sumerr, flag)
it = np.broadcast(x, y, a, b, theta, r, rin, rout, seg_id, sum, sumerr, flag)
while np.PyArray_MultiIter_NOTDONE(it):
status = sep_sum_ellipse(
Expand Down Expand Up @@ -1404,7 +1409,14 @@ def sum_ellipse(np.ndarray data not None, x, y, a, b, theta, r=1.0,
subpix, 0, &bkgflux, &bkgfluxerr, &bkgarea, &bkgflag)
_assert_ok(status)

if area1 > 0:
if not bkgarea > 0:
raise ValueError(
"The background annulus does not contain any valid pixels, "
"for the object at index "
f"{np.PyArray_MultiIter_INDEX(it)}."
)

if (area1 > 0):
flux1 -= bkgflux / bkgarea * area1
bkgfluxerr = bkgfluxerr / bkgarea * area1
fluxerr1 = sqrt(fluxerr1*fluxerr1 + bkgfluxerr*bkgfluxerr)
Expand Down