/usr/lib/python2.7/dist-packages/photutils/isophote/tests/test_regression.py is in python-photutils 0.4-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 | # Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Despite being cast as a unit test, this code implements regression
testing of the Ellipse algorithm, against results obtained by the
stsdas$analysis/isophote task 'ellipse'.
The stsdas task was run on test images and results were stored in
tables. The code here runs the Ellipse algorithm on the same images,
producing a list of Isophote instances. The contents of this list then
get compared with the contents of the corresponding table.
Some quantities are compared in assert statements. These were designed
to be executed only when the synth_highsnr.fits image is used as input.
That way, we are mainly checking numerical differences that originate in
the algorithms themselves, and not caused by noise. The quantities
compared this way are:
- mean intensity: less than 1% diff. for sma > 3 pixels, 5% otherwise
- ellipticity: less than 1% diff. for sma > 3 pixels, 20% otherwise
- position angle: less than 1 deg. diff. for sma > 3 pixels, 20 deg.
otherwise
- X and Y position: less than 0.2 pixel diff.
For the M51 image we have mostly good agreement with the SPP code in
most of the parameters (mean isophotal intensity agrees within a
fraction of 1% mostly), but every now and then the ellipticity and
position angle of the semi-major axis may differ by a large amount from
what the SPP code measures. The code also stops prematurely wrt the
larger sma values measured by the SPP code. This is caused by a
difference in the way the gradient relative error is measured in each
case, and suggests that the SPP code may have a bug.
The not-so-good behavior observed in the case of the M51 image is to be
expected though. This image is exactly the type of galaxy image for
which the algorithm *wasn't* designed for. It has an almost negligible
smooth ellipsoidal component, and a lot of lumpy spiral structure that
causes the radial gradient computation to go berserk. On top of that,
the ellipticity is small (roundish isophotes) throughout the image,
causing large relative errors and instability in the fitting algorithm.
For now, we can only check the bilinear integration mode. The mean and
median modes cannot be checked since the original 'ellipse' task has a
bug that causes the creation of erroneous output tables. A partial
comparison could be made if we write new code that reads the standard
output of 'ellipse' instead, captured from screen, and use it as
reference for the regression.
"""
from __future__ import (absolute_import, division, print_function,
unicode_literals)
import math
import numpy as np
import os.path as op
import pytest
from astropy.io import fits
from astropy.table import Table
from astropy.tests.helper import remote_data
from ..ellipse import Ellipse
from ..integrator import BILINEAR
from ...datasets import get_path
try:
import scipy # noqa
HAS_SCIPY = True
except ImportError:
HAS_SCIPY = False
@remote_data
@pytest.mark.skipif('not HAS_SCIPY')
# @pytest.mark.parametrize('name', ['M51', 'synth', 'synth_lowsnr',
# 'synth_highsnr'])
@pytest.mark.parametrize('name', ['synth_highsnr'])
def test_regression(name, integrmode=BILINEAR, verbose=False):
"""
NOTE: The original code in SPP won't create the right table for the MEAN
integration moder, so use the screen output at synth_table_mean.txt to
compare results visually with synth_table_mean.fits.
"""
filename = '{0}_table.fits'.format(name)
path = op.join(op.dirname(op.abspath(__file__)), 'data', filename)
table = Table.read(path)
nrows = len(table['SMA'])
path = get_path('isophote/{0}.fits'.format(name),
location='photutils-datasets', cache=True)
hdu = fits.open(path)
data = hdu[0].data
hdu.close()
ellipse = Ellipse(data)
isophote_list = ellipse.fit_image()
# isophote_list = ellipse.fit_image(sclip=2., nclip=3)
fmt = ("%5.2f %6.1f %8.3f %8.3f %8.3f %9.5f %6.2f "
"%6.2f %6.2f %5.2f %4d %3d %3d %2d")
for row in range(nrows):
try:
iso = isophote_list[row]
except IndexError:
# skip non-existent rows in isophote list, if that's the case.
break
# data from Isophote
sma_i = iso.sample.geometry.sma
intens_i = iso.intens
int_err_i = iso.int_err if iso.int_err else 0.
pix_stddev_i = iso.pix_stddev if iso.pix_stddev else 0.
rms_i = iso.rms if iso.rms else 0.
ellip_i = iso.sample.geometry.eps if iso.sample.geometry.eps else 0.
pa_i = iso.sample.geometry.pa if iso.sample.geometry.pa else 0.
x0_i = iso.sample.geometry.x0
y0_i = iso.sample.geometry.y0
rerr_i = (iso.sample.gradient_relative_error
if iso.sample.gradient_relative_error else 0.)
ndata_i = iso.ndata
nflag_i = iso.nflag
niter_i = iso.niter
stop_i = iso.stop_code
# convert to old code reference system
pa_i = (pa_i - np.pi/2) / np.pi * 180.
x0_i += 1
y0_i += 1
# ref data from table
sma_t = table['SMA'][row]
intens_t = table['INTENS'][row]
int_err_t = table['INT_ERR'][row]
pix_stddev_t = table['PIX_VAR'][row]
rms_t = table['RMS'][row]
ellip_t = table['ELLIP'][row]
pa_t = table['PA'][row]
x0_t = table['X0'][row]
y0_t = table['Y0'][row]
rerr_t = table['GRAD_R_ERR'][row]
ndata_t = table['NDATA'][row]
nflag_t = table['NFLAG'][row]
niter_t = table['NITER'][row] if table['NITER'][row] else 0
stop_t = table['STOP'][row] if table['STOP'][row] else -1
# relative differences
sma_d = (sma_i - sma_t) / sma_t * 100. if sma_t > 0. else 0.
intens_d = (intens_i - intens_t) / intens_t * 100.
int_err_d = ((int_err_i - int_err_t) / int_err_t * 100.
if int_err_t > 0. else 0.)
pix_stddev_d = ((pix_stddev_i - pix_stddev_t) / pix_stddev_t * 100.
if pix_stddev_t > 0. else 0.)
rms_d = (rms_i - rms_t) / rms_t * 100. if rms_t > 0. else 0.
ellip_d = (ellip_i - ellip_t) / ellip_t * 100.
pa_d = pa_i - pa_t # diff in angle is absolute
x0_d = x0_i - x0_t # diff in position is absolute
y0_d = y0_i - y0_t
rerr_d = rerr_i - rerr_t # diff in relative error is absolute
ndata_d = (ndata_i - ndata_t) / ndata_t * 100.
nflag_d = 0
niter_d = 0
stop_d = 0 if stop_i == stop_t else -1
if verbose:
print("* data " + fmt % (sma_i, intens_i, int_err_i, pix_stddev_i,
rms_i, ellip_i, pa_i, x0_i, y0_i, rerr_i,
ndata_i, nflag_i, niter_i, stop_i))
print(" ref " + fmt % (sma_t, intens_t, int_err_t, pix_stddev_t,
rms_t, ellip_t, pa_t, x0_t, y0_t, rerr_t,
ndata_t, nflag_t, niter_t, stop_t))
print(" diff " + fmt % (sma_d, intens_d, int_err_d, pix_stddev_d,
rms_d, ellip_d, pa_d, x0_d, y0_d, rerr_d,
ndata_d, nflag_d, niter_d, stop_d))
print()
if name == "synth_highsnr" and integrmode == BILINEAR:
assert abs(x0_d) <= 0.21
assert abs(y0_d) <= 0.21
if sma_i > 3.:
assert abs(intens_d) <= 1.
else:
assert abs(intens_d) <= 5.
if not math.isnan(ellip_d):
if sma_i > 3.:
assert abs(ellip_d) <= 1. # 1%
else:
assert abs(ellip_d) <= 20. # 20%
if not math.isnan(pa_d):
if sma_i > 3.:
assert abs(pa_d) <= 1. # 1 deg.
else:
assert abs(pa_d) <= 20. # 20 deg.
|