forked from sep-developers/sep
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsep.pyx
2136 lines (1769 loc) · 76.4 KB
/
sep.pyx
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
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# This wrapper licensed under an MIT license.
"""
Source Extraction and Photometry
This module is a wrapper of the SEP C library.
"""
import numpy as np
cimport numpy as np
from libc cimport limits
from libc.math cimport sqrt
cimport cython
from cpython.mem cimport PyMem_Malloc, PyMem_Free
from cpython.version cimport PY_MAJOR_VERSION
np.import_array() # To access the numpy C-API.
__version__ = "1.2.1"
# -----------------------------------------------------------------------------
# Definitions from the SEP C library
# macro definitions from sep.h
DEF SEP_TBYTE = 11
DEF SEP_TINT = 31
DEF SEP_TFLOAT = 42
DEF SEP_TDOUBLE = 82
# input flag values (C macros)
DEF SEP_NOISE_NONE = 0
DEF SEP_NOISE_STDDEV = 1
DEF SEP_NOISE_VAR = 2
# filter types for sep_extract
DEF SEP_FILTER_CONV = 0
DEF SEP_FILTER_MATCHED = 1
# Threshold types
DEF SEP_THRESH_REL = 0
DEF SEP_THRESH_ABS = 1
# input flags for aperture photometry
DEF SEP_MASK_IGNORE = 0x0004
# Output flag values accessible from python
OBJ_MERGED = np.short(0x0001)
OBJ_TRUNC = np.short(0x0002)
OBJ_DOVERFLOW = np.short(0x0004)
OBJ_SINGU = np.short(0x0008)
APER_TRUNC = np.short(0x0010)
APER_HASMASKED = np.short(0x0020)
APER_ALLMASKED = np.short(0x0040)
APER_NONPOSITIVE = np.short(0x0080)
# macro defintion from sepcore.h
# This is not part of the SEP API, but we pull it out because want to
# explicitly detect memory errors so that we can raise MemoryError().
DEF MEMORY_ALLOC_ERROR = 1
# header definitions
cdef extern from "sep.h":
ctypedef struct sep_image:
const void *data
const void *noise
const void *mask
const void *segmap
int dtype
int ndtype
int mdtype
int sdtype
int w
int h
double noiseval
short noise_type
double gain
double maskthresh
ctypedef struct sep_bkg:
int w
int h
float globalback
float globalrms
ctypedef struct sep_catalog:
int nobj
float *thresh
int *npix
int *tnpix
int *xmin
int *xmax
int *ymin
int *ymax
double *x
double *y
double *x2
double *y2
double *xy
double *errx2
double *erry2
double *errxy
float *a
float *b
float *theta
float *cxx
float *cyy
float *cxy
float *cflux
float *flux
float *cpeak
float *peak
int *xcpeak
int *ycpeak
int *xpeak
int *ypeak
short *flag
int **pix
int *objectspix
int sep_background(const sep_image *im,
int bw, int bh,
int fw, int fh,
double fthresh,
sep_bkg **bkg)
float sep_bkg_global(const sep_bkg *bkg)
float sep_bkg_globalrms(const sep_bkg *bkg)
int sep_bkg_array(const sep_bkg *bkg, void *arr, int dtype)
int sep_bkg_rmsarray(const sep_bkg *bkg, void *arr, int dtype)
int sep_bkg_subarray(const sep_bkg *bkg, void *arr, int dtype)
void sep_bkg_free(sep_bkg *bkg)
int sep_extract(const sep_image *image,
float thresh,
int thresh_type,
int minarea,
float *conv,
int convw, int convh,
int filter_type,
int deblend_nthresh,
double deblend_cont,
int clean_flag,
double clean_param,
sep_catalog **catalog)
void sep_catalog_free(sep_catalog *catalog)
int sep_sum_circle(const sep_image *image,
double x, double y, double r,
int id, int subpix, short inflags,
double *sum, double *sumerr, double *area, short *flag)
int sep_sum_circann(const sep_image *image,
double x, double y, double rin, double rout,
int id, int subpix, short inflags,
double *sum, double *sumerr, double *area, short *flag)
int sep_sum_ellipse(const sep_image *image,
double x, double y, double a, double b, double theta,
double r, int id, int subpix, short inflags,
double *sum, double *sumerr, double *area,
short *flag)
int sep_sum_ellipann(const sep_image *image,
double x, double y, double a, double b,
double theta, double rin, double rout,
int id, int subpix,
short inflags,
double *sum, double *sumerr, double *area,
short *flag)
int sep_flux_radius(const sep_image *image,
double x, double y, double rmax, int id, int subpix,
short inflag,
double *fluxtot, double *fluxfrac, int n,
double *r, short *flag)
int sep_kron_radius(const sep_image *image,
double x, double y, double cxx, double cyy,
double cxy, double r, int id,
double *kronrad, short *flag)
int sep_windowed(const sep_image *image,
double x, double y, double sig, int subpix, short inflag,
double *xout, double *yout, int *niter, short *flag)
int sep_ellipse_axes(double cxx, double cyy, double cxy,
double *a, double *b, double *theta)
void sep_ellipse_coeffs(double a, double b, double theta,
double *cxx, double *cyy, double *cxy)
void sep_set_ellipse(unsigned char *arr, int w, int h,
double x, double y,
double cxx, double cyy, double cxy, double r,
unsigned char val)
void sep_set_extract_pixstack(size_t val)
size_t sep_get_extract_pixstack()
void sep_set_sub_object_limit(int val)
int sep_get_sub_object_limit()
void sep_get_errmsg(int status, char *errtext)
void sep_get_errdetail(char *errtext)
# -----------------------------------------------------------------------------
# Utility functions
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. "
"Only native byte order arrays are supported. "
"To change the byte order of the array 'data', do "
"'data = data.byteswap().newbyteorder()'".format(dtype))
t = dtype.type
if t is np.single:
return SEP_TFLOAT
elif t is np.bool_ or t is np.ubyte:
return SEP_TBYTE
elif dtype == np.double:
return SEP_TDOUBLE
elif dtype == np.intc:
return SEP_TINT
raise ValueError('input array dtype not supported: {0}'.format(dtype))
cdef int _check_array_get_dims(np.ndarray arr, int *w, int *h) except -1:
"""Check some things about an array and return dimensions"""
# Raise an informative message if array is not C-contiguous
if not arr.flags["C_CONTIGUOUS"]:
raise ValueError("array is not C-contiguous")
# Check that there are exactly 2 dimensions
if arr.ndim != 2:
raise ValueError("array must be 2-d")
# ensure that arr dimensions are not too large for C ints.
if arr.shape[0] <= <Py_ssize_t> limits.INT_MAX:
h[0] = arr.shape[0]
else:
raise ValueError("array height ({0:d}) greater than INT_MAX ({1:d})"
.format(arr.shape[0], limits.INT_MAX))
if arr.shape[1] <= <Py_ssize_t> limits.INT_MAX:
w[0] = arr.shape[1]
else:
raise ValueError("array width ({0:d}) greater than INT_MAX ({1:d})"
.format(arr.shape[1], limits.INT_MAX))
return 0
cdef int _assert_ok(int status) except -1:
"""Get the SEP error message corresponding to status code"""
cdef char *errmsg
cdef char *errdetail
if status == 0:
return 0
# First check if we have an out-of-memory error, so we don't try to
# allocate more memory to hold the error message.
if status == MEMORY_ALLOC_ERROR:
raise MemoryError
# Otherwise, get error message.
errmsg = <char *>PyMem_Malloc(61 * sizeof(char))
sep_get_errmsg(status, errmsg)
pyerrmsg = <bytes> errmsg
PyMem_Free(errmsg)
# Get error detail.
errdetail = <char *>PyMem_Malloc(512 * sizeof(char))
sep_get_errdetail(errdetail)
pyerrdetail = <bytes> errdetail
PyMem_Free(errdetail)
# If error detail is present, append it to the message.
if pyerrdetail != b"":
pyerrmsg = pyerrmsg + b": " + pyerrdetail
# Convert string to unicode if on python 3
if PY_MAJOR_VERSION == 3:
msg = pyerrmsg.decode()
else:
msg = pyerrmsg
raise Exception(msg)
cdef int _parse_arrays(np.ndarray data, err, var, mask, segmap,
sep_image *im) except -1:
"""Helper function for functions accepting data, error, mask & segmap arrays.
Fills in an sep_image struct."""
cdef int ew, eh, mw, mh, sw, sh
cdef np.uint8_t[:,:] buf, ebuf, mbuf, sbuf
# Clear im fields we might not touch (everything besides data, dtype, w, h)
im.noise = NULL
im.mask = NULL
im.segmap = NULL
im.ndtype = 0
im.mdtype = 0
im.noiseval = 0.0
im.noise_type = SEP_NOISE_NONE
im.gain = 0.0
im.maskthresh = 0.0
# Get main image info
_check_array_get_dims(data, &(im.w), &(im.h))
im.dtype = _get_sep_dtype(data.dtype)
buf = data.view(dtype=np.uint8)
im.data = <void*>&buf[0, 0]
# Check if noise is error or variance.
noise = None # will point to either error or variance.
if err is not None:
if var is not None:
raise ValueError("Cannot specify both err and var")
noise = err
im.noise_type = SEP_NOISE_STDDEV
elif var is not None:
noise = var
im.noise_type = SEP_NOISE_VAR
# parse noise
if noise is None:
im.noise = NULL
im.noise_type = SEP_NOISE_NONE
im.noiseval = 0.0
elif isinstance(noise, np.ndarray):
if noise.ndim == 0:
im.noise = NULL
im.noiseval = noise
elif noise.ndim == 2:
_check_array_get_dims(noise, &ew, &eh)
if ew != im.w or eh != im.h:
raise ValueError("size of error/variance array must match"
" data")
im.ndtype = _get_sep_dtype(noise.dtype)
ebuf = noise.view(dtype=np.uint8)
im.noise = <void*>&ebuf[0, 0]
else:
raise ValueError("error/variance array must be 0-d or 2-d")
else:
im.noise = NULL
im.noiseval = noise
# Optional input: mask
if mask is None:
im.mask = NULL
else:
_check_array_get_dims(mask, &mw, &mh)
if mw != im.w or mh != im.h:
raise ValueError("size of mask array must match data")
im.mdtype = _get_sep_dtype(mask.dtype)
mbuf = mask.view(dtype=np.uint8)
im.mask = <void*>&mbuf[0, 0]
# Optional input: segmap
if segmap is None:
im.segmap = NULL
else:
_check_array_get_dims(segmap, &sw, &sh)
if sw != im.w or sh != im.h:
raise ValueError("size of segmap array must match data")
im.sdtype = _get_sep_dtype(segmap.dtype)
sbuf = segmap.view(dtype=np.uint8)
im.segmap = <void*>&sbuf[0, 0]
# -----------------------------------------------------------------------------
# Background Estimation
cdef class Background:
"""
Background(data, mask=None, maskthresh=0.0, bw=64, bh=64,
fw=3, fh=3, fthresh=0.0)
Representation of spatially variable image background and noise.
Parameters
----------
data : 2-d `~numpy.ndarray`
Data array.
mask : 2-d `~numpy.ndarray`, optional
Mask array, optional
maskthresh : float, optional
Mask threshold. This is the inclusive upper limit on the mask value
in order for the corresponding pixel to be unmasked. For boolean
arrays, False and True are interpreted as 0 and 1, respectively.
Thus, given a threshold of zero, True corresponds to masked and
False corresponds to unmasked.
bw, bh : int, optional
Size of background boxes in pixels. Default is 64.
fw, fh : int, optional
Filter width and height in boxes. Default is 3.
fthresh : float, optional
Filter threshold. Default is 0.0.
"""
cdef sep_bkg *ptr # pointer to C struct
cdef np.dtype orig_dtype # dtype code of original image
@cython.boundscheck(False)
@cython.wraparound(False)
def __cinit__(self, np.ndarray data not None, np.ndarray mask=None,
float maskthresh=0.0, int bw=64, int bh=64,
int fw=3, int fh=3, float fthresh=0.0):
cdef int status
cdef sep_image im
_parse_arrays(data, None, None, mask, None, &im)
im.maskthresh = maskthresh
status = sep_background(&im, bw, bh, fw, fh, fthresh, &self.ptr)
_assert_ok(status)
self.orig_dtype = data.dtype
# Note: all initialization work is done in __cinit__. This is just here
# for the docstring.
def __init__(self, np.ndarray data not None, np.ndarray mask=None,
float maskthresh=0.0, int bw=64, int bh=64,
int fw=3, int fh=3, float fthresh=0.0):
"""Background(data, mask=None, maskthresh=0.0, bw=64, bh=64,
fw=3, fh=3, fthresh=0.0)"""
pass
property globalback:
"""Global background level."""
def __get__(self):
return sep_bkg_global(self.ptr)
property globalrms:
"""Global background RMS."""
def __get__(self):
return sep_bkg_globalrms(self.ptr)
def back(self, dtype=None):
"""back(dtype=None)
Create an array of the background.
Parameters
----------
dtype : `~numpy.dtype`, optional
Data type of output array. Default is the dtype of the original
data.
Returns
-------
back : `~numpy.ndarray`
Array with same dimensions as original data.
"""
cdef int sep_dtype
cdef np.uint8_t[:, :] buf
if dtype is None:
dtype = self.orig_dtype
else:
dtype = np.dtype(dtype)
sep_dtype = _get_sep_dtype(dtype)
result = np.empty((self.ptr.h, self.ptr.w), dtype=dtype)
buf = result.view(dtype=np.uint8)
status = sep_bkg_array(self.ptr, &buf[0, 0], sep_dtype)
_assert_ok(status)
return result
def rms(self, dtype=None):
"""rms(dtype=None)
Create an array of the background rms.
Parameters
----------
dtype : `~numpy.dtype`, optional
Data type of output array. Default is the dtype of the original
data.
Returns
-------
rms : `~numpy.ndarray`
Array with same dimensions as original data.
"""
cdef int sep_dtype
cdef np.uint8_t[:, :] buf
if dtype is None:
dtype = self.orig_dtype
else:
dtype = np.dtype(dtype)
sep_dtype = _get_sep_dtype(dtype)
result = np.empty((self.ptr.h, self.ptr.w), dtype=dtype)
buf = result.view(dtype=np.uint8)
status = sep_bkg_rmsarray(self.ptr, &buf[0, 0], sep_dtype)
_assert_ok(status)
return result
def subfrom(self, np.ndarray data not None):
"""subfrom(data)
Subtract the background from an existing array.
Like ``data = data - bkg``, but avoids making a copy of the data.
Parameters
----------
data : `~numpy.ndarray`
Input array, which will be updated in-place. Shape must match
that of the original image used to measure the background.
"""
cdef int w, h, status, sep_dtype
cdef np.uint8_t[:, :] buf
assert self.ptr is not NULL
_check_array_get_dims(data, &w, &h)
sep_dtype = _get_sep_dtype(data.dtype)
buf = data.view(dtype=np.uint8)
# ensure dimensions match original image
if (w != self.ptr.w or h != self.ptr.h):
raise ValueError("Data dimensions do not match background "
"dimensions")
status = sep_bkg_subarray(self.ptr, &buf[0, 0], sep_dtype)
_assert_ok(status)
def __array__(self, dtype=None):
return self.back(dtype=dtype)
def __rsub__(self, np.ndarray data not None):
data = np.copy(data)
self.subfrom(data)
return data
def __dealloc__(self):
if self.ptr is not NULL:
sep_bkg_free(self.ptr)
# -----------------------------------------------------------------------------
# Source Extraction
# This needs to match the result from extract
cdef packed struct Object:
np.float64_t thresh
np.int_t npix
np.int_t tnpix
np.int_t xmin
np.int_t xmax
np.int_t ymin
np.int_t ymax
np.float64_t x
np.float64_t y
np.float64_t x2
np.float64_t y2
np.float64_t xy
np.float64_t a
np.float64_t b
np.float64_t theta
np.float64_t cxx
np.float64_t cyy
np.float64_t cxy
np.float64_t cflux
np.float64_t flux
np.float64_t cpeak
np.float64_t peak
np.float64_t errx2
np.float64_t erry2
np.float64_t errxy
np.int_t xcpeak
np.int_t ycpeak
np.int_t xpeak
np.int_t ypeak
np.int_t flag
default_kernel = np.array([[1.0, 2.0, 1.0],
[2.0, 4.0, 2.0],
[1.0, 2.0, 1.0]], dtype=np.float32)
def extract(np.ndarray data not None, float thresh, err=None, var=None,
gain=None, np.ndarray mask=None, double maskthresh=0.0,
int minarea=5,
np.ndarray filter_kernel=default_kernel, filter_type='matched',
int deblend_nthresh=32, double deblend_cont=0.005,
bint clean=True, double clean_param=1.0,
segmentation_map=False):
"""extract(data, thresh, err=None, mask=None, minarea=5,
filter_kernel=default_kernel, filter_type='matched',
deblend_nthresh=32, deblend_cont=0.005, clean=True,
clean_param=1.0, segmentation_map=False)
Extract sources from an image.
Parameters
----------
data : `~numpy.ndarray`
Data array (2-d).
thresh : float
Threshold pixel value for detection. If an ``err`` or ``var`` array
is not given, this is interpreted as an absolute threshold. If ``err``
or ``var`` is given, this is interpreted as a relative threshold: the
absolute threshold at pixel (j, i) will be ``thresh * err[j, i]`` or
``thresh * sqrt(var[j, i])``.
err, var : float or `~numpy.ndarray`, optional
Error *or* variance (specify at most one). This can be used to
specify a pixel-by-pixel detection threshold; see "thresh" argument.
gain : float, optional
Conversion factor between data array units and poisson counts. This
does not affect detection; it is used only in calculating Poisson
noise contribution to uncertainty parameters such as ``errx2``. If
not given, no Poisson noise will be added.
mask : `~numpy.ndarray`, optional
Mask array. ``True`` values, or numeric values greater than
``maskthresh``, are considered masked. Masking a pixel is equivalent
to setting data to zero and noise (if present) to infinity.
maskthresh : float, optional
Threshold for a pixel to be masked. Default is ``0.0``.
minarea : int, optional
Minimum number of pixels required for an object. Default is 5.
filter_kernel : `~numpy.ndarray` or None, optional
Filter kernel used for on-the-fly filtering (used to
enhance detection). Default is a 3x3 array:
[[1,2,1], [2,4,2], [1,2,1]]. Set to ``None`` to skip
convolution.
filter_type : {'matched', 'conv'}, optional
Filter treatment. This affects filtering behavior when a noise
array is supplied. ``'matched'`` (default) accounts for
pixel-to-pixel noise in the filter kernel. ``'conv'`` is
simple convolution of the data array, ignoring pixel-to-pixel
noise across the kernel. ``'matched'`` should yield better
detection of faint sources in areas of rapidly varying noise
(such as found in coadded images made from semi-overlapping
exposures). The two options are equivalent when noise is
constant.
deblend_nthresh : int, optional
Number of thresholds used for object deblending. Default is 32.
deblend_cont : float, optional
Minimum contrast ratio used for object deblending. Default is 0.005.
To entirely disable deblending, set to 1.0.
clean : bool, optional
Perform cleaning? Default is True.
clean_param : float, optional
Cleaning parameter (see SExtractor manual). Default is 1.0.
segmentation_map : bool, optional
If True, also return a "segmentation map" giving the member
pixels of each object. Default is False.
Returns
-------
objects : `~numpy.ndarray`
Extracted object parameters (structured array). Available fields are:
* ``thresh`` (float) Threshold at object location.
* ``npix`` (int) Number of pixels belonging to the object.
* ``tnpix`` (int) Number of pixels above threshold (unconvolved data).
* ``xmin``, ``xmax`` (int) Minimum, maximum x coordinates of pixels.
* ``ymin``, ``ymax`` (int) Minimum, maximum y coordinates of pixels.
* ``x``, ``y`` (float) object barycenter (first moments).
* ``x2``, ``y2``, ``xy`` (float) Second moments.
* ``errx2``, ``erry2``, ``errxy`` (float) Second moment errors.
Note that these will be zero if error is not given.
* ``a``, ``b``, ``theta`` (float) Ellipse parameters, scaled as
described by Section 8.4.2 in "The Source Extractor Guide" or
Section 10.1.5-6 of v2.13 of SExtractor's User Manual.
* ``cxx``, ``cyy``, ``cxy`` (float) Alternative ellipse parameters.
* ``cflux`` (float) Sum of member pixels in convolved data.
* ``flux`` (float) Sum of member pixels in unconvolved data.
* ``cpeak`` (float) Peak value in convolved data.
* ``peak`` (float) Peak value in unconvolved data.
* ``xcpeak``, ``ycpeak`` (int) Coordinate of convolved peak pixel.
* ``xpeak``, ``ypeak`` (int) Coordinate of unconvolved peak pixel.
* ``flag`` (int) Extraction flags.
segmap : `~numpy.ndarray`, optional
Array of integers with same shape as data. Pixels not belonging to
any object have value 0. All pixels belonging to the ``i``-th object
(e.g., ``objects[i]``) have value ``i+1``. Only returned if
``segmentation_map=True``.
"""
cdef int kernelw, kernelh, status, i, j
cdef int filter_typecode, thresh_type
cdef sep_catalog *catalog = NULL
cdef np.ndarray[Object] result
cdef float[:, :] kernelflt
cdef float *kernelptr
cdef np.int32_t[:, :] segmap_buf
cdef np.int32_t *segmap_ptr
cdef int *objpix
cdef sep_image im
# parse arrays
_parse_arrays(data, err, var, mask, None, &im)
im.maskthresh = maskthresh
if gain is not None:
im.gain = gain
# Parse filter input
if filter_kernel is None:
kernelptr = NULL
kernelw = 0
kernelh = 0
else:
kernelflt = filter_kernel.astype(np.float32)
kernelptr = &kernelflt[0, 0]
kernelw = kernelflt.shape[1]
kernelh = kernelflt.shape[0]
if filter_type == 'matched':
filter_typecode = SEP_FILTER_MATCHED
elif filter_type == 'conv':
filter_typecode = SEP_FILTER_CONV
else:
raise ValueError("unknown filter_type: {!r}".format(filter_type))
# If image has error info, the threshold is relative, otherwise
# it is absolute.
if im.noise_type == SEP_NOISE_NONE:
thresh_type = SEP_THRESH_ABS
else:
thresh_type = SEP_THRESH_REL
status = sep_extract(&im,
thresh, thresh_type, minarea,
kernelptr, kernelw, kernelh, filter_typecode,
deblend_nthresh, deblend_cont, clean, clean_param,
&catalog)
_assert_ok(status)
# Allocate result record array and fill it
result = np.empty(catalog.nobj,
dtype=np.dtype([('thresh', np.float64),
('npix', np.int_),
('tnpix', np.int_),
('xmin', np.int_),
('xmax', np.int_),
('ymin', np.int_),
('ymax', np.int_),
('x', np.float64),
('y', np.float64),
('x2', np.float64),
('y2', np.float64),
('xy', np.float64),
('errx2', np.float64),
('erry2', np.float64),
('errxy', np.float64),
('a', np.float64),
('b', np.float64),
('theta', np.float64),
('cxx', np.float64),
('cyy', np.float64),
('cxy', np.float64),
('cflux', np.float64),
('flux', np.float64),
('cpeak', np.float64),
('peak', np.float64),
('xcpeak', np.int_),
('ycpeak', np.int_),
('xpeak', np.int_),
('ypeak', np.int_),
('flag', np.int_)]))
for i in range(catalog.nobj):
result['thresh'][i] = catalog.thresh[i]
result['npix'][i] = catalog.npix[i]
result['tnpix'][i] = catalog.tnpix[i]
result['xmin'][i] = catalog.xmin[i]
result['xmax'][i] = catalog.xmax[i]
result['ymin'][i] = catalog.ymin[i]
result['ymax'][i] = catalog.ymax[i]
result['x'][i] = catalog.x[i]
result['y'][i] = catalog.y[i]
result['x2'][i] = catalog.x2[i]
result['y2'][i] = catalog.y2[i]
result['xy'][i] = catalog.xy[i]
result['errx2'][i] = catalog.errx2[i]
result['erry2'][i] = catalog.erry2[i]
result['errxy'][i] = catalog.errxy[i]
result['a'][i] = catalog.a[i]
result['b'][i] = catalog.b[i]
result['theta'][i] = catalog.theta[i]
result['cxx'][i] = catalog.cxx[i]
result['cyy'][i] = catalog.cyy[i]
result['cxy'][i] = catalog.cxy[i]
result['cflux'][i] = catalog.cflux[i]
result['flux'][i] = catalog.flux[i]
result['cpeak'][i] = catalog.cpeak[i]
result['peak'][i] = catalog.peak[i]
result['xcpeak'][i] = catalog.xcpeak[i]
result['ycpeak'][i] = catalog.ycpeak[i]
result['xpeak'][i] = catalog.xpeak[i]
result['ypeak'][i] = catalog.ypeak[i]
result['flag'][i] = catalog.flag[i]
# construct a segmentation map, if it was requested.
if segmentation_map:
# Note: We have to write out `(data.shape[0], data.shape[1])` because
# because Cython turns `data.shape` later into an int pointer when
# the function argument is typed as np.ndarray.
segmap = np.zeros((data.shape[0], data.shape[1]), dtype=np.int32)
segmap_buf = segmap
segmap_ptr = &segmap_buf[0, 0]
for i in range(catalog.nobj):
objpix = catalog.pix[i]
for j in range(catalog.npix[i]):
segmap_ptr[objpix[j]] = i + 1
# Free the C catalog
sep_catalog_free(catalog)
if segmentation_map:
return result, segmap
else:
return result
# -----------------------------------------------------------------------------
# Aperture Photometry
@cython.boundscheck(False)
@cython.wraparound(False)
def sum_circle(np.ndarray data not None, x, y, r,
var=None, err=None, gain=None, np.ndarray mask=None,
double maskthresh=0.0,
seg_id=None, np.ndarray segmap=None,
bkgann=None, int subpix=5):
"""sum_circle(data, x, y, r, err=None, var=None, mask=None, maskthresh=0.0,
segmap=None, seg_id=None,
bkgann=None, gain=None, subpix=5)
Sum data in circular aperture(s).
Parameters
----------
data : `~numpy.ndarray`
2-d array to be summed.
x, y, r : array_like
Center coordinates and radius (radii) of aperture(s).
``x`` corresponds to the second ("fast") axis of the input array
and ``y`` corresponds to the first ("slow") axis.
``x, y = (0.0, 0.0)`` corresponds to the center of the first
element of the array. These inputs obey numpy broadcasting rules.
err, var : float or `~numpy.ndarray`
Error *or* variance (specify at most one).
mask : `~numpy.ndarray`, optional
Mask array. If supplied, a given pixel is masked if its value
is greater than ``maskthresh``.
maskthresh : float, optional
Threshold for a pixel to be masked. Default is ``0.0``.
segmap : `~numpy.ndarray`, optional
Segmentation image with dimensions of ``data`` and dtype ``np.int32``.
This is an optional input and corresponds to the segmentation map
output by `~sep.extract`.
seg_id : array_like, optional
Array of segmentation ids used to mask additional pixels in the image.
Dimensions correspond to the dimensions of ``x`` and ``y``. The
behavior differs depending on whether ``seg_id`` is negative or
positive. If ``seg_id`` is positive, all pixels belonging to other
objects are masked. (Pixel ``j, i`` is masked if ``seg[j, i] != seg_id
and seg[j, i] != 0``). If ``seg_id`` is negative, all pixels other
than those belonging to the object of interest are masked. (Pixel ``j,
i`` is masked if ``seg[j, i] != -seg_id``). NB: must be included if
``segmap` is provided.
bkgann : tuple, optional
Length 2 tuple giving the inner and outer radius of a
"background annulus". If supplied, the background is estimated
by averaging unmasked pixels in this annulus. If supplied, the inner
and outer radii obey numpy broadcasting rules along with ``x``,
``y`` and ``r``.
gain : float, optional
Conversion factor between data array units and poisson counts,
used in calculating poisson noise in aperture sum. If ``None``
(default), do not add poisson noise.
subpix : int, optional
Subpixel sampling factor. If 0, exact overlap is calculated.
Default is 5.
Returns
-------
sum : `~numpy.ndarray`
The sum of the data array within the aperture.
sumerr : `~numpy.ndarray`
Error on the sum.
flags : `~numpy.ndarray`
Integer giving flags. (0 if no flags set.)
"""
cdef double flux1, fluxerr1, area1,
cdef double bkgflux, bkgfluxerr, bkgarea
cdef short flag1, bkgflag
cdef size_t i
cdef int status
cdef np.broadcast it
cdef sep_image im
# Test for map without seg_id. Nothing happens if seg_id supplied but
# without segmap.
if (segmap is not None) and (seg_id is None):
raise ValueError('`segmap` supplied but not `seg_id`.')
_parse_arrays(data, err, var, mask, segmap, &im)
im.maskthresh = maskthresh
if gain is not None:
im.gain = gain
# Require that inputs are float64 arrays. This has to be done because we
# are using a broadcasting iterator below, where we need to know the type
# in advance. There are other ways to do this, e.g., using NpyIter_Multi
# in the numpy C-API. However, the best way to use this from cython
# is not clear to me at this time.
#
# docs.scipy.org/doc/numpy/reference/c-api.iterator.html#NpyIter_MultiNew
dt = np.dtype(np.double)
dint = np.dtype(np.int32)
x = np.require(x, dtype=dt)
y = np.require(y, dtype=dt)
r = np.require(r, dtype=dt)
# Segmentation image and ids with same dimensions as x, y, etc.
if seg_id is not None:
seg_id = np.require(seg_id, dtype=dint)
if seg_id.shape != x.shape:
raise ValueError('Shapes of `x` and `seg_id` do not match')
else:
seg_id = np.zeros(len(x), dtype=dint)
if bkgann is None:
# allocate ouput arrays
shape = np.broadcast(x, y, r).shape
sum = np.empty(shape, dt)
sumerr = np.empty(shape, dt)
flag = np.empty(shape, np.short)
it = np.broadcast(x, y, r, seg_id, sum, sumerr, flag)
while np.PyArray_MultiIter_NOTDONE(it):
status = sep_sum_circle(
&im,
(<double*>np.PyArray_MultiIter_DATA(it, 0))[0],
(<double*>np.PyArray_MultiIter_DATA(it, 1))[0],
(<double*>np.PyArray_MultiIter_DATA(it, 2))[0],
(<int*>np.PyArray_MultiIter_DATA(it, 3))[0],
subpix, 0,
<double*>np.PyArray_MultiIter_DATA(it, 4),
<double*>np.PyArray_MultiIter_DATA(it, 5),
&area1,
<short*>np.PyArray_MultiIter_DATA(it, 6))
_assert_ok(status)
# Advance the iterator
np.PyArray_MultiIter_NEXT(it)
return sum, sumerr, flag
else:
rin, rout = bkgann
# Require float arrays (see note above)
rin = np.require(rin, dtype=dt)
rout = np.require(rout, dtype=dt)
# allocate ouput arrays
shape = np.broadcast(x, y, r, rin, rout).shape
sum = np.empty(shape, dt)
sumerr = np.empty(shape, dt)
flag = np.empty(shape, np.short)
it = np.broadcast(x, y, r, rin, rout, seg_id, sum, sumerr, flag)
while np.PyArray_MultiIter_NOTDONE(it):
status = sep_sum_circle(
&im,
(<double*>np.PyArray_MultiIter_DATA(it, 0))[0],
(<double*>np.PyArray_MultiIter_DATA(it, 1))[0],
(<double*>np.PyArray_MultiIter_DATA(it, 2))[0],
(<int*>np.PyArray_MultiIter_DATA(it, 5))[0],
subpix, 0, &flux1, &fluxerr1, &area1, &flag1)
_assert_ok(status)