Skip to content

Instantly share code, notes, and snippets.

@emmeowzing
Last active August 28, 2019 02:10
Show Gist options
  • Save emmeowzing/8a7de265bb7ebe2384dca3bd4350668b to your computer and use it in GitHub Desktop.
Save emmeowzing/8a7de265bb7ebe2384dca3bd4350668b to your computer and use it in GitHub Desktop.
Generate Gaussian and Moffat filters for SExtractor
# Default configuration file for SExtractor 2.19.5
# EB 2013-12-14
#
#-------------------------------- Catalog ------------------------------------
PARAMETERS_NAME output.param # name of the file containing catalog contents
#------------------------------- Extraction ----------------------------------
DETECT_TYPE CCD # CCD (linear) or PHOTO (with gamma correction)
FILTER Y # apply filter for detection (Y or N)?
FILTER_NAME filters/filter_6x6.conv # name of the file containing the filter
SEEING_FWHM 1.987 # stellar FWHM in arcsec
# I arrived at this value after running my dist-fit
# program and multiplying the output FWHM value by
# the size of an individual pixel in arcseconds
# (in this case 0.74, since the camera was binned
# 1x1)
STARNNW_NAME default.nnw # Neural-Network_Weight table filename
#------------------------------- Photometry ----------------------------------
PHOT_APERTURES 10 # MAG_APER aperture diameter(s) in pixels
PHOT_AUTOPARAMS 2.5, 3.5 # MAG_AUTO parameters: <Kron_fact>,<min_radius>
PHOT_PETROPARAMS 2.0, 3.5 # MAG_PETRO parameters: <Petrosian_fact>,
# <min_radius>
SATUR_LEVEL 60000.0 # level (in ADUs) at which arises saturation
SATUR_KEY SATURATE # keyword for saturation level (in ADUs)
MAG_ZEROPOINT 15.0 # magnitude zero-point corrective factor
GAIN 138.0 # detector gain in e-/ADU * exptime
GAIN_KEY GAIN # keyword for detector gain in e-/ADU
PIXEL_SCALE 0.74 # size of pixel in arcsec (0=use FITS WCS info)
BACK_SIZE 64 # Background mesh: <size> or <width>,<height>
BACK_FILTERSIZE 3 # Background filter: <size> or <width>,<height>
BACKPHOTO_TYPE GLOBAL # can be GLOBAL or LOCAL
#--------------------- Memory (change with caution!) -------------------------
MEMORY_OBJSTACK 3000 # number of objects in stack
MEMORY_PIXSTACK 1000000 # number of pixels in stack
MEMORY_BUFSIZE 1024 # number of lines in buffer
#----------------------------- Miscellaneous ---------------------------------
VERBOSE_TYPE QUIET # can be QUIET, NORMAL or FULL
#! /usr/bin/env python3.6
# -*- coding: utf-8 -*-
""" Command line script to fabricate SExtractor filters (either Moffat or Gaussian) """
from typing import Tuple, List
from types import FunctionType as Function
from functools import wraps
from scipy.integrate import dblquad
from math import pi
import numpy as np
import argparse
import sys
def scale(f: Function) -> Function:
""" Decorator to conserve information """
@wraps(f)
def _new_f(*args, **kwargs) -> np.ndarray:
arr = f(*args, **kwargs)
# Scale array by volume (in this case \sum_{i,j}k_{i,j} in the discrete case
# is the same as a double integral to find volume on
# \iint_{[-3 * sigma, 3 * sigma]} in the continuous case/
return arr / np.sum(arr)
return _new_f
class Kernel:
""" Generate different kernels for testing """
def __init__(self, dims: Tuple[int, int]) -> None:
self.dims = dims
@scale
def gaussian(self, sigma: float =1.0) -> np.ndarray:
kernel = np.empty(self.dims)
unit_square = (-0.5, 0.5, lambda y: -0.5, lambda y: 0.5)
x_shift = 0 if self.dims[0] % 2 else 0.5
y_shift = 0 if self.dims[1] % 2 else 0.5
for i in range(self.dims[0]):
for j in range(self.dims[1]):
# integrate on a unit square centered at the origin as the
# function moves about it in discrete unit steps
res = dblquad(
lambda x, y: 1 / (2 * pi * sigma ** 2) * np.exp(
- ((x + i - self.dims[0] // 2 + x_shift) / sigma) ** 2
- ((y + j - self.dims[1] // 2 + y_shift) / sigma) ** 2),
*unit_square
)[0]
kernel[i][j] = res
return kernel
@scale
def moffat(self, alpha: float =3.0, beta: float =2.5) -> np.ndarray:
""" The Gaussian is a limiting case of this kernel as \beta -> \infty """
kernel = np.zeros(self.dims)
unit_square = (-0.5, 0.5, lambda y: -0.5, lambda y: 0.5)
x_shift = 0 if self.dims[0] % 2 else 0.5
y_shift = 0 if self.dims[1] % 2 else 0.5
for i in range(self.dims[0]):
for j in range(self.dims[1]):
res = dblquad(
lambda x, y: (beta - 1) / (pi * alpha ** 2) * \
(1 + ((x + i - self.dims[0] // 2 + x_shift) ** 2
+ (y + j - self.dims[1] // 2 + y_shift) ** 2)
/ (alpha ** 2)) ** (-beta),
*unit_square
)[0]
kernel[i][j] = res
return kernel
def _write(arr: np.ndarray, output_filename: str, params: List[float]) -> None:
# Use NumPy to save this array to file with proper formatting
np.savetxt(output_filename, arr, fmt='%.6f',
header=f"{'Gaussian ' if len(params) == 1 else 'Moffat '} \
convolution filter for params {', '.join(map(str, params))}"
)
if __name__ == '__main__':
class ErrParser(argparse.ArgumentParser):
def error_msg(self, message: str) -> None:
sys.stderr.write(f'Error: {message}\n')
self.print_help()
sys.exit(1)
parser = ErrParser(description=__doc__)
## create a Gaussian or Moffat filter
group = parser.add_mutually_exclusive_group(required=True)
group.add_argument(
'--gaussian',
action='store_true',
help='Create a Gaussian filter'
)
group.add_argument(
'--moffat',
action='store_true',
help='Create a Moffat filter'
)
parser.add_argument(
'-i',
'--input_config',
type=str,
help='Input config file with parameters; '
'if none is given, default parameters are assigned'
)
parser.add_argument(
'-s',
'--size',
type=int,
default=5,
help='Size of the output filter'
)
parser.add_argument(
'-o',
'--output_file',
type=str,
default='filters/filter_5x5.conv',
help='Output filter filename'
)
args = parser.parse_args()
# handle various configurations
if args.size < 3:
parser.error(f'Expected requested output filter size to be >2, received {args.size}')
kernel = Kernel((args.size, args.size))
with open(args.input_config) as file:
try:
params = list(map(float, file))
except ValueError:
print('error: Must supply at least one value in input configuration file')
sys.exit(1)
kern = np.zeros((args.size,)*2)
if args.gaussian:
try:
kern = kernel.gaussian(*params)
except TypeError:
parser.error('Expected 1 argument to output Gaussian convolution filter,'
f'received {len(params)}')
else:
try:
kern = kernel.moffat(*params)
except TypeError:
parser.error('Expected 2 arguments to output Moffat convolution filter,'
f'received {len(params)}')
_write(kern, args.output_file, params)
#np.set_printoptions(precision=2, suppress=True)
#print(kern)
#! /bin/bash
# Run SExtractor on an image
# Print how to use this script
usage()
{
echo "Usage: $0 -I images" 1>&2
exit 1
}
# Check for an input list of files using getopts
while getopts "I:" opt
do
case $opt in
I)
IMAGES=${OPTARG}
# Check that input images exist
for image in $IMAGES
do
if ! [ -e $image ]
then
echo "** Image $image does not exist, exiting" 1>&2
exit 1
fi
done
;;
*)
usage
;;
esac
done
# Make this argument a required one
if [ -z "${IMAGES}" ]
then
usage
fi
# VARIABLES
SIGMA=2.686 # Determine this with dist-fitting
GAUSSIAN_CONFIG=gaussian_config.conf # Config containing gaussian params
DIM=6 # Dimension square Gaussian filter
FILTER=filters/filter_${DIM}x${DIM}.conv # location of Gaussian filter
touch $FILTER
if ! [ -e $GAUSSIAN_CONFIG ]
then
echo $SIGMA > $GAUSSIAN_CONFIG
fi
# Generate the Gaussian convolution filter using psf.py
echo "** Generating Gaussian *.conv filter"
./psf.py --gaussian -i gaussian_config.conf -s $DIM -o $FILTER
# Insert `CONV NORM` at the beginning of the file, indicating that
echo -e "CONV NORM\n$(cat $FILTER)" > $FILTER
# Now set up SExtractor for star extraction
# If this file doesn't exist we can create it quite easily
create_output_param()
{
# We need this info for differential photometry;
# see `sextractr -dp` for a list of output vals
echo "NUMBER" > output.param
echo "X_IMAGE" >> output.param
echo "Y_IMAGE" >> output.param
echo "MAG_BEST" >> output.param
echo "MAGERR_AUTO" >> output.param
}
if ! [ -e output.param ]
then
create_output_param
fi
# Run SE on the list of images
echo "** Running SExtractor"
for image in $IMAGES
do
echo "** Extracting sources from $image"
output="${image%.*}_data"
# make a directory for storing data if it doesn't exist
if ! [ -e $output ]
then
mkdir $output
fi
# run SE on this image
sextractor $image -c config.sex
# sort the output file by magnitude (i.e. column 4)
sort -k 4 -n "test.cat" > "new_test.cat"
rm test.cat
mv new_test.cat $output/sources.dat
done
exit 0
# 1 NUMBER Running object number
# 2 X_IMAGE Object position along x [pixel]
# 3 Y_IMAGE Object position along y [pixel]
# 4 MAG_BEST Best of MAG_AUTO and MAG_ISOCOR [mag]
# 5 MAGERR_AUTO RMS error for AUTO magnitude [mag]
37 377.0922 259.7026 1.8010 0.0009
53 550.5630 179.5938 1.8824 0.0010
24 748.1370 332.1368 3.3383 0.0036
54 557.9897 178.7481 3.6416 0.0044
7 461.4656 509.2312 3.8440 0.0036
14 108.1221 432.5323 4.4696 0.0092
47 267.4219 234.8703 5.2240 0.0167
40 479.7549 262.6016 5.5575 0.0243
6 60.4464 463.3252 5.8106 0.0287
17 329.1845 414.0296 5.8514 0.0310
71 537.5078 97.7215 5.9777 0.0335
76 759.6041 78.0504 6.0854 0.0387
32 510.6487 295.6030 6.0964 0.0374
23 350.6108 347.7325 6.1269 0.0451
62 185.1494 166.4438 6.1352 0.0373
59 266.1955 171.5730 6.2263 0.0387
41 408.6657 258.4180 6.2314 0.0433
31 52.5987 298.8428 6.3077 0.0628
33 675.3002 280.5554 6.3093 0.0483
4 425.0814 6.8543 6.4334 0.0578
9 579.7260 461.7381 6.5671 0.0714
39 77.9612 261.9085 6.5958 0.0581
30 329.1181 307.3435 6.7043 0.0588
38 459.3319 266.8521 6.7471 0.0843
45 166.5703 240.5071 6.7993 0.0603
18 347.4777 410.6228 6.8219 0.0775
63 561.8380 148.9742 6.8405 0.0833
1 316.2875 14.5624 6.9140 0.0816
78 470.3698 65.0490 7.0081 0.0855
87 338.9956 492.4057 7.0119 0.1042
25 313.3618 335.0482 7.0148 0.0950
19 736.5134 401.4340 7.0159 0.0711
13 347.1508 441.4143 7.0244 0.0873
11 604.3749 446.2296 7.0405 0.1126
20 655.7521 369.9928 7.0767 0.0759
27 468.5163 322.0381 7.0943 0.0920
80 716.0372 46.1597 7.1686 0.1178
86 673.8177 498.7320 7.1948 0.1170
58 27.1819 172.7696 7.2227 0.0837
56 599.2001 180.5753 7.2528 0.1322
16 437.9594 430.6689 7.3005 0.1446
73 421.7666 90.7760 7.3057 0.1272
89 749.7230 474.2617 7.3593 0.1264
22 75.7799 349.8125 7.4120 0.0968
69 354.5099 111.3286 7.4591 0.0938
72 292.2825 91.9467 7.4754 0.1188
29 698.7561 312.7304 7.4776 0.1504
82 321.0089 38.1161 7.4879 0.1028
84 338.0222 22.8832 7.5475 0.1405
85 475.0169 16.7784 7.5628 0.1167
67 10.2820 126.3758 7.5888 0.1468
83 645.1759 23.5960 7.6148 0.1504
2 229.0993 14.5086 7.6217 0.1024
46 725.1545 238.6853 7.6264 0.1685
35 343.4024 274.0893 7.6869 0.1597
81 573.0203 43.4791 7.7111 0.1572
61 141.2247 171.2217 7.7219 0.1301
10 42.2475 453.4137 7.7223 0.1423
8 154.4733 509.1797 7.7300 0.0954
43 330.4054 247.3587 7.7340 0.1805
79 485.5110 47.7709 7.7471 0.1919
74 80.5979 90.9078 7.7515 0.1509
15 641.7231 433.0136 7.7901 0.1756
55 300.8346 182.8186 7.7988 0.1140
3 368.9584 11.5983 7.7998 0.1327
64 396.6917 145.1164 7.8248 0.1921
52 390.4153 201.0126 7.8273 0.1461
21 242.9762 364.9773 7.8641 0.1937
26 583.9088 335.1248 7.8712 0.1904
50 744.4726 221.2071 7.8810 0.1430
77 265.3370 73.7371 7.9028 0.1948
75 686.2951 89.9559 7.9212 0.1957
49 704.1933 222.4802 7.9439 0.1802
68 674.6395 124.3786 7.9949 0.1946
12 607.7104 442.0107 8.0067 0.2211
5 31.7227 4.8798 8.0134 0.1616
65 642.5117 134.5186 8.0294 0.1256
48 486.9579 234.3194 8.1299 0.2002
34 103.2579 282.2443 8.1541 0.1993
51 531.2029 213.1771 8.1912 0.2024
60 76.9658 173.5108 8.2067 0.1845
57 82.0089 178.2706 8.2395 0.1879
36 764.6843 272.4473 8.2462 0.1416
66 603.9733 127.6699 8.2485 0.1894
88 522.1910 480.1825 8.2504 0.2074
44 499.9854 243.0162 8.2668 0.2040
42 217.9463 249.5039 8.4811 0.1466
28 436.1946 318.4135 8.5137 0.2589
70 431.1633 103.8057 8.5711 0.2550
# 1 NUMBER Running object number
# 2 X_IMAGE Object position along x [pixel]
# 3 Y_IMAGE Object position along y [pixel]
# 4 MAG_BEST Best of MAG_AUTO and MAG_ISOCOR [mag]
# 5 MAGERR_AUTO RMS error for AUTO magnitude [mag]
37 379.7175 259.7540 1.8115 0.0009
53 553.2457 179.6278 1.8841 0.0010
25 750.7803 332.1775 3.3440 0.0038
54 560.6536 178.8930 3.7433 0.0042
7 463.9523 509.1723 3.8460 0.0040
14 110.7914 432.5323 4.4958 0.0099
46 270.0834 234.9817 5.2549 0.0176
38 482.3987 262.6208 5.5301 0.0239
6 63.0111 463.4528 5.7769 0.0296
17 332.0390 414.0395 5.8012 0.0329
70 540.1830 97.7062 5.9643 0.0332
23 353.2594 347.9436 6.1096 0.0450
75 762.1459 78.1858 6.1290 0.0365
33 513.3340 295.5502 6.1605 0.0400
61 187.7349 166.5515 6.1902 0.0399
58 268.9426 171.5294 6.2069 0.0395
41 411.2591 258.3527 6.2421 0.0500
34 677.7612 280.7247 6.2691 0.0506
9 582.2932 461.8884 6.3281 0.0737
32 55.3242 298.8999 6.4262 0.0514
4 427.9588 7.0338 6.5570 0.0592
39 80.6646 262.1248 6.6859 0.0569
12 606.9709 446.1689 6.7449 0.0929
31 331.5746 307.2919 6.7935 0.0687
44 169.0456 240.5943 6.8036 0.0648
36 461.8581 266.8176 6.8659 0.0857
19 739.1589 401.4252 6.8702 0.0684
63 564.5526 149.0272 6.8932 0.0862
13 349.7422 441.6076 6.9422 0.0676
18 350.1137 411.0287 6.9764 0.0698
78 472.6908 65.0210 7.0121 0.0891
86 676.7299 498.6837 7.0449 0.1224
26 315.6629 334.9966 7.0564 0.0978
1 318.8411 14.4852 7.0646 0.0731
21 658.2527 369.9692 7.0701 0.0880
50 520.6581 217.5814 7.0826 0.0743
28 471.1284 322.1149 7.1000 0.0905
56 602.2910 181.0930 7.1169 0.1232
72 423.9444 90.8253 7.2298 0.1203
80 718.4812 46.3767 7.2715 0.1267
24 78.2921 349.7142 7.3188 0.1182
16 440.2726 430.5268 7.3588 0.1473
59 30.2861 172.6294 7.3917 0.0903
89 752.8344 474.1380 7.4500 0.1257
60 144.2459 171.2392 7.4754 0.1432
76 311.9620 73.2659 7.4768 0.1272
85 477.8297 17.0741 7.4862 0.1453
87 341.8714 492.5576 7.4917 0.0924
11 44.8296 453.3908 7.5255 0.1373
73 83.1142 91.2604 7.5371 0.1033
71 294.8715 92.1287 7.5577 0.1222
81 576.0319 43.4994 7.5641 0.1126
8 157.3616 509.2314 7.5973 0.1215
3 371.6983 11.6235 7.6069 0.1488
69 357.1154 111.2008 7.6134 0.1377
83 647.4932 23.4982 7.6248 0.1808
22 245.7028 364.7695 7.6505 0.1741
2 231.6921 14.5391 7.6610 0.1010
30 701.2806 312.8606 7.6816 0.1556
52 393.1328 200.9527 7.7200 0.1475
35 346.1532 273.9568 7.7528 0.1791
82 323.5526 38.2539 7.7824 0.1631
57 80.0830 173.8168 7.7980 0.1754
77 267.8461 73.8593 7.7988 0.1928
43 333.0092 246.9980 7.8661 0.1662
68 677.5017 124.2247 7.9140 0.1963
84 340.9916 23.1921 7.9156 0.1258
40 454.7679 262.4005 7.9237 0.1845
67 12.9153 126.1656 7.9327 0.1712
49 746.8461 221.1191 7.9429 0.2454
65 642.4686 134.4666 7.9507 0.1221
55 303.4991 183.4715 7.9993 0.1865
5 34.3525 5.0673 8.0662 0.1936
64 399.5290 145.4818 8.0934 0.2187
88 525.2077 480.6071 8.1438 0.2130
42 220.6625 249.9111 8.1525 0.1470
48 706.6909 222.4098 8.1725 0.2028
74 690.1339 89.8340 8.1858 0.2179
15 643.0591 432.5051 8.2107 0.2335
47 489.9438 234.5377 8.2388 0.1971
29 437.9819 318.4843 8.2475 0.2537
66 606.9835 127.7188 8.2508 0.2405
51 534.2907 213.0174 8.2911 0.2761
45 728.0050 238.5046 8.3560 0.2196
27 586.0206 334.9911 8.4196 0.2088
62 704.8127 164.3910 8.4197 0.2198
79 488.8218 47.8416 8.4284 0.2590
20 184.7257 392.5778 8.4453 0.2050
10 471.9806 461.5148 8.5170 0.2519
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment