/usr/lib/python3/dist-packages/photutils/isophote/model.py is in python3-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 | # Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import (absolute_import, division, print_function,
unicode_literals)
import numpy as np
from .geometry import EllipseGeometry
__all__ = ['build_ellipse_model']
def build_ellipse_model(shape, isolist, fill=0., high_harmonics=False):
"""
Build an elliptical model 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
the isophote "pixel". The algorithm takes care of this partial
pixel coverage by keeping track of how much intensity was added to
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.
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`
The image with the model galaxy.
"""
from scipy.interpolate import LSQUnivariateSpline
# the target grid is spaced in 0.1 pixel intervals so as
# to ensure no gaps will result on the output array.
finely_spaced_sma = np.arange(isolist[0].sma, isolist[-1].sma, 0.1)
# interpolate ellipse parameters
# End points must be discarded, but how many?
# This seems to work so far
nodes = isolist.sma[2:-2]
intens_array = LSQUnivariateSpline(
isolist.sma, isolist.intens, nodes)(finely_spaced_sma)
eps_array = LSQUnivariateSpline(
isolist.sma, isolist.eps, nodes)(finely_spaced_sma)
pa_array = LSQUnivariateSpline(
isolist.sma, isolist.pa, nodes)(finely_spaced_sma)
x0_array = LSQUnivariateSpline(
isolist.sma, isolist.x0, nodes)(finely_spaced_sma)
y0_array = LSQUnivariateSpline(
isolist.sma, isolist.y0, nodes)(finely_spaced_sma)
grad_array = LSQUnivariateSpline(
isolist.sma, isolist.grad, nodes)(finely_spaced_sma)
a3_array = LSQUnivariateSpline(
isolist.sma, isolist.a3, nodes)(finely_spaced_sma)
b3_array = LSQUnivariateSpline(
isolist.sma, isolist.b3, nodes)(finely_spaced_sma)
a4_array = LSQUnivariateSpline(
isolist.sma, isolist.a4, nodes)(finely_spaced_sma)
b4_array = LSQUnivariateSpline(
isolist.sma, isolist.b4, nodes)(finely_spaced_sma)
# Return deviations from ellipticity to their original amplitude meaning
a3_array = -a3_array * grad_array * finely_spaced_sma
b3_array = -b3_array * grad_array * finely_spaced_sma
a4_array = -a4_array * grad_array * finely_spaced_sma
b4_array = -b4_array * grad_array * finely_spaced_sma
# 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 outside image boundaries, ignore.
if (i > 0 and i < shape[0]-1 and j > 0 and j < shape[1] - 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 = geometry.radius(phi)
# zero weight values must be set to 1.
weight[np.where(weight <= 0.)] = 1.
# normalize
result /= weight
# fill value
result[np.where(result == 0.)] = fill
return result
|