Skip to content

Commit 6093f49

Browse files
timstaleytimstaley
committed
Initial implementation of moments-estimation and source-fitting.
2 parents 9eeda90 + 3ccd31b commit 6093f49

20 files changed

+2234
-174
lines changed

‎docs/source/notebooks.rst

+2-1
Original file line numberDiff line numberDiff line change
@@ -7,4 +7,5 @@ Jupyter Notebook Examples
77
notebooks/convolutional_kernel_funcs
88
notebooks/anti_aliasing_properties
99
notebooks/simulate_and_image
10-
notebooks/sourcefinding
10+
notebooks/sourcefinding
11+
notebooks/gaussian_parameterizations

‎docs/source/notebooks/gaussian_parameterizations.ipynb

+276
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,276 @@
1+
{
2+
"cells": [
3+
{
4+
"cell_type": "markdown",
5+
"metadata": {},
6+
"source": [
7+
"# Re-parameterising 2d Gaussian Profiles"
8+
]
9+
},
10+
{
11+
"cell_type": "markdown",
12+
"metadata": {},
13+
"source": [
14+
"In the general case of a multivariate Gaussian, the PDF is as follows:"
15+
]
16+
},
17+
{
18+
"cell_type": "markdown",
19+
"metadata": {},
20+
"source": [
21+
"$$\n",
22+
"f_{\\mathbf X}(x_1,\\ldots,x_k) = \\frac{\\exp\\left(-\\frac 1 2 ({\\mathbf x}-{\\boldsymbol\\mu})^\\mathrm{T}{\\boldsymbol\\Sigma}^{-1}({\\mathbf x}-{\\boldsymbol\\mu})\\right)}{\\sqrt{(2\\pi)^k|\\boldsymbol\\Sigma|}}\n",
23+
"$$"
24+
]
25+
},
26+
{
27+
"cell_type": "markdown",
28+
"metadata": {},
29+
"source": [
30+
"For the 2d case, setting"
31+
]
32+
},
33+
{
34+
"cell_type": "markdown",
35+
"metadata": {},
36+
"source": [
37+
"$$\n",
38+
"\\boldsymbol\\mu = \\begin{pmatrix} \\mu_X \\\\ \\mu_Y \\end{pmatrix}, \\quad\n",
39+
" \\boldsymbol\\Sigma = \\begin{pmatrix} \\sigma_X^2 & \\rho \\sigma_X \\sigma_Y \\\\\n",
40+
" \\rho \\sigma_X \\sigma_Y & \\sigma_Y^2 \\end{pmatrix}\n",
41+
"$$"
42+
]
43+
},
44+
{
45+
"cell_type": "markdown",
46+
"metadata": {},
47+
"source": [
48+
"so that \n",
49+
"$$\n",
50+
" |\\boldsymbol\\Sigma| = \\sigma_X^2 \\sigma_Y^2(1 - \\rho^2)\n",
51+
"$$\n",
52+
"\n",
53+
"and \n",
54+
"\n",
55+
"$$\n",
56+
"\\boldsymbol\\Sigma^{-1} = \\frac{1}{|\\boldsymbol\\Sigma|} \\begin{pmatrix} \n",
57+
" \\sigma_Y^2 & -\\rho \\sigma_X \\sigma_Y \\\\\n",
58+
" -\\rho \\sigma_X \\sigma_Y & \\sigma_X^2 \n",
59+
" \\end{pmatrix}\n",
60+
"$$"
61+
]
62+
},
63+
{
64+
"cell_type": "markdown",
65+
"metadata": {},
66+
"source": [
67+
"Then substituting and simplifying gives [Eq 1]:\n",
68+
"$$\n",
69+
"\\begin{align}\n",
70+
" f(x,y) &=\n",
71+
" \\frac{1}{2 \\pi \\sigma_X \\sigma_Y \\sqrt{1-\\rho^2}}\n",
72+
" \\exp\\left(\n",
73+
" -\\frac{1}{2(1-\\rho^2)}\\left[\n",
74+
" \\frac{(x-\\mu_X)^2}{\\sigma_X^2} +\n",
75+
" \\frac{(y-\\mu_Y)^2}{\\sigma_Y^2} -\n",
76+
" \\frac{2\\rho(x-\\mu_X)(y-\\mu_Y)}{\\sigma_X \\sigma_Y}\n",
77+
" \\right]\n",
78+
" \\right)\\\\ \n",
79+
"\\end{align}\n",
80+
"$$"
81+
]
82+
},
83+
{
84+
"cell_type": "markdown",
85+
"metadata": {},
86+
"source": [
87+
"For a source profile, the normalising constant factor is simply replaced by $A$, the peak amplitude.\n",
88+
"If there is no correlation ($\\rho =0$, which can alternatively be viewed as the major/minor profile axes being aligned with x/y) then this simplifies to the 'separable' 2D-Gaussian:\n",
89+
"$$\n",
90+
"\\begin{align}\n",
91+
" f(x,y) &=\n",
92+
" A\n",
93+
" \\exp\\left(\n",
94+
" -\\frac{1}{2}\\left[\n",
95+
" \\frac{(x-\\mu_X)^2}{\\sigma_X^2} +\n",
96+
" \\frac{(y-\\mu_Y)^2}{\\sigma_Y^2}\n",
97+
" \\right]\n",
98+
" \\right)\\\\\n",
99+
"\\end{align}\n",
100+
"$$"
101+
]
102+
},
103+
{
104+
"cell_type": "markdown",
105+
"metadata": {},
106+
"source": [
107+
"Alternatively we can use the following parameterisation, where $\\theta$ represents the rotation of the major-axis of the Gaussian *counterclockwise* relative to the x-positive axis, and $\\sigma_{maj}$, $\\sigma_{min}$ represent the equivalent standard deviations in x and y if the Gaussian were *not rotated* - (i.e. the semimajor and semiminor axes) [Eq 2]:"
108+
]
109+
},
110+
{
111+
"cell_type": "markdown",
112+
"metadata": {},
113+
"source": [
114+
"$$\n",
115+
"f(x,y) = A \\exp\\left(- \\left(a(x - \\mu_X)^2 + b(x-\\mu_X)(y-\\mu_Y) + c(y-\\mu_Y)^2 \\right)\\right)\n",
116+
"$$"
117+
]
118+
},
119+
{
120+
"cell_type": "markdown",
121+
"metadata": {},
122+
"source": [
123+
"Where\n",
124+
"$$\n",
125+
"a = \\frac{1}{2}\\left(\\frac{\\cos^2\\theta}{\\sigma_{maj}^2} + \\frac{\\sin^2\\theta}{\\sigma_{min}^2}\\right) \\\\\n",
126+
"b = \\frac{1}{2}\\left(\\frac{\\sin2\\theta}{\\sigma_{maj}^2} + \\frac{\\sin2\\theta}{\\sigma_{min}^2}\\right) \\\\\n",
127+
"c = \\frac{1}{2}\\left(\\frac{\\sin^2\\theta}{\\sigma_{maj}^2} + \\frac{\\cos^2\\theta}{\\sigma_{min}^2}\\right) \\\\\n",
128+
"$$"
129+
]
130+
},
131+
{
132+
"cell_type": "markdown",
133+
"metadata": {},
134+
"source": [
135+
"(This is the opposite rotation to that shown on Wikipedia, and matches the Astropy implementation.)\n",
136+
"\n",
137+
"This is convenient when representing source profiles - we can immediately read off the relative size of the semi-major and semi-minor axes, interpret the source rotation angle, etc. Note, $\\sigma_{maj}$ as defined here is the standard deviation equivalent length - to get the FWHM of the source profile you must multiply this by a factor of $2\\sqrt{2\\ln{2}}$ (see e.g. [Mathworld](http://mathworld.wolfram.com/GaussianFunction.html) for a derivation). (This matches the implementations in Astropy and S-Extractor, but PySE quotes the HWHM, i.e. larger by a factor $\\sqrt{2\\ln{2}}$.)\n",
138+
"\n",
139+
"However, the rotation-angle parameterisation makes comparison of source fits trickier. For $\\mu$ and $\\sigma$ values we can use standard 'are these floating point values close' routines, but $\\theta$ presents problems. For a start, even if we confine it to the range $(-\\pi,\\pi)$, we can still get two approximate fits with reported rotation angles just inside each of those bounds which actually represent very nearly identical fits. Second, when $\\sigma_x=\\sigma_y$, the rotation angle $\\theta$ is effectively degenerate (since the profile is circular), so we can ignore it in comparisons. How do we smoothly transition between ignoring it completely and taking it into account for elongated profiles, when comparing noisy fits?\n",
140+
"\n",
141+
"It seems a better approach would be to compare covariance matrices in the reference x-y frame. That way we effectively replace $\\theta$ with the correlation $\\rho$, which varies smoothly in the range $(-1,1)$, and simply tends to zero as a the Gaussian circularizes.\n",
142+
"\n",
143+
"To do so, we take advantage of the fact that the covariance matrix for the rotated Gaussian is simply:"
144+
]
145+
},
146+
{
147+
"cell_type": "markdown",
148+
"metadata": {},
149+
"source": [
150+
"$$\n",
151+
"\\Sigma_{rot} = \\begin{pmatrix} \\sigma_{maj}^2 & 0 \\\\\n",
152+
" 0 & \\sigma_{min}^2 \\end{pmatrix}\n",
153+
"$$\n",
154+
"\n",
155+
"in the rotated frame. To transform to the reference frame we define the rotation matrix:\n",
156+
"\n",
157+
"$$\n",
158+
"R = \\begin{pmatrix} cos(\\theta) & sin(\\theta) \\\\\n",
159+
" -sin(\\theta) & \\cos(\\theta) \\end{pmatrix}\n",
160+
"$$\n",
161+
"\n",
162+
"and apply a change of basis transform:\n",
163+
"\n",
164+
"$$\n",
165+
"\\Sigma_{ref} = R\\,\\Sigma_{rot}\\,R^{-1}\n",
166+
"$$\n",
167+
"\n",
168+
"Then we compare reference-frame covariance matrices (or specifically, correlation coeffcients) for the source fits in question."
169+
]
170+
},
171+
{
172+
"cell_type": "code",
173+
"execution_count": null,
174+
"metadata": {},
175+
"outputs": [],
176+
"source": [
177+
"%load_ext autoreload\n",
178+
"%autoreload 2"
179+
]
180+
},
181+
{
182+
"cell_type": "code",
183+
"execution_count": null,
184+
"metadata": {},
185+
"outputs": [],
186+
"source": [
187+
"from fastimgproto.sourcefind.fit import Gaussian2dParams\n",
188+
"import numpy as np\n",
189+
"import astropy.units as u\n",
190+
"import attr"
191+
]
192+
},
193+
{
194+
"cell_type": "code",
195+
"execution_count": null,
196+
"metadata": {},
197+
"outputs": [],
198+
"source": [
199+
"from astropy.coordinates import Angle, SkyCoord"
200+
]
201+
},
202+
{
203+
"cell_type": "code",
204+
"execution_count": null,
205+
"metadata": {},
206+
"outputs": [],
207+
"source": [
208+
"g1 = Gaussian2dParams(x_centre=0,\n",
209+
" y_centre=0,\n",
210+
" amplitude=1,\n",
211+
" semimajor=1.5,\n",
212+
" semiminor=1,\n",
213+
" theta=np.pi/4,\n",
214+
" )\n",
215+
"\n",
216+
"g1"
217+
]
218+
},
219+
{
220+
"cell_type": "code",
221+
"execution_count": null,
222+
"metadata": {},
223+
"outputs": [],
224+
"source": [
225+
"g1.covariance"
226+
]
227+
},
228+
{
229+
"cell_type": "code",
230+
"execution_count": null,
231+
"metadata": {},
232+
"outputs": [],
233+
"source": [
234+
"g1.correlation"
235+
]
236+
},
237+
{
238+
"cell_type": "code",
239+
"execution_count": null,
240+
"metadata": {},
241+
"outputs": [],
242+
"source": [
243+
"g1 = attr.evolve(g1, theta = -np.pi/4)\n",
244+
"g1.correlation"
245+
]
246+
},
247+
{
248+
"cell_type": "code",
249+
"execution_count": null,
250+
"metadata": {},
251+
"outputs": [],
252+
"source": [
253+
"g1 = attr.evolve(g1, theta = np.pi/4, semimajor=1000)\n",
254+
"g1.correlation"
255+
]
256+
},
257+
{
258+
"cell_type": "code",
259+
"execution_count": null,
260+
"metadata": {},
261+
"outputs": [],
262+
"source": [
263+
"g1 = attr.evolve(g1, semimajor=1.0001)\n",
264+
"g1.correlation"
265+
]
266+
}
267+
],
268+
"metadata": {
269+
"language_info": {
270+
"name": "python",
271+
"pygments_lexer": "ipython3"
272+
}
273+
},
274+
"nbformat": 4,
275+
"nbformat_minor": 1
276+
}

‎docs/source/notebooks/simulate_and_image.ipynb

+2-2
Original file line numberDiff line numberDiff line change
@@ -241,8 +241,8 @@
241241
"ax1.set_xlim(axlims)\n",
242242
"ax1.set_ylim(axlims)\n",
243243
"for src in sfimage.islands:\n",
244-
" ax1.axvline(src.xbar, ls=':')\n",
245-
" ax1.axhline(src.ybar, ls=':')"
244+
" ax1.axvline(src.params.moments_fit.x_centre, ls=':')\n",
245+
" ax1.axhline(src.params.moments_fit.y_centre, ls=':')"
246246
]
247247
}
248248
],

‎docs/source/notebooks/sourcefinding.ipynb

+7-5
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,7 @@
7676
"srcs.append(fixture.gaussian_point_source(x_centre=64.12, y_centre=48.88, amplitude=1.0))\n",
7777
"srcs.append(fixture.gaussian_point_source(x_centre=128.43, y_centre=94.5, amplitude=1.5))\n",
7878
"for s in srcs:\n",
79-
" img_data += fixture.evaluate_model_on_pixel_grid(img_data.shape, s)"
79+
" fixture.add_gaussian2d_to_image(s, img_data)"
8080
]
8181
},
8282
{
@@ -185,9 +185,10 @@
185185
"outputs": [],
186186
"source": [
187187
"plt.imshow(island.data)\n",
188-
"plt.xlim(island.extremum_x_idx-10,island.extremum_x_idx+10,)\n",
189-
"plt.ylim(island.extremum_y_idx-10,island.extremum_y_idx+10,)\n",
190-
"plt.scatter(island.xbar,island.ybar, marker='*', s=200, c='y',)"
188+
"plt.xlim(island.extremum.index.x-10,island.extremum.index.x+10,)\n",
189+
"plt.ylim(island.extremum.index.y-10,island.extremum.index.y+10,)\n",
190+
"moments_fit = island.params.moments_fit\n",
191+
"plt.scatter(moments_fit.x_centre,moments_fit.y_centre, marker='*', s=200, c='y',)"
191192
]
192193
},
193194
{
@@ -201,7 +202,8 @@
201202
"print()\n",
202203
"print(\"Island barycentres:\")\n",
203204
"for i in sfimage.islands:\n",
204-
" print(i.xbar, i.ybar)"
205+
" moments = i.params.moments_fit\n",
206+
" print(moments.x_centre, moments.y_centre)"
205207
]
206208
}
207209
],

0 commit comments

Comments
 (0)