File indexing completed on 2024-04-28 15:10:12
0001 /* 0002 This file is part of SEP 0003 0004 SPDX-FileCopyrightText: 1993-2011 Emmanuel Bertin -- IAP /CNRS/UPMC 0005 SPDX-FileCopyrightText: 2014 SEP developers 0006 0007 SPDX-License-Identifier: LGPL-3.0-or-later 0008 */ 0009 0010 /* datatype codes */ 0011 #define SEP_TBYTE 11 /* 8-bit unsigned byte */ 0012 #define SEP_TINT 31 /* native int type */ 0013 #define SEP_TFLOAT 42 0014 #define SEP_TDOUBLE 82 0015 0016 /* object & aperture flags */ 0017 #define SEP_OBJ_MERGED 0x0001 /* object is result of deblending */ 0018 #define SEP_OBJ_TRUNC 0x0002 /* object truncated at image boundary */ 0019 #define SEP_OBJ_DOVERFLOW 0x0004 /* not currently used, but could be */ 0020 #define SEP_OBJ_SINGU 0x0008 /* x,y fully correlated */ 0021 #define SEP_APER_TRUNC 0x0010 0022 #define SEP_APER_HASMASKED 0x0020 0023 #define SEP_APER_ALLMASKED 0x0040 0024 #define SEP_APER_NONPOSITIVE 0x0080 0025 0026 /* noise_type values in sep_image */ 0027 #define SEP_NOISE_NONE 0 0028 #define SEP_NOISE_STDDEV 1 0029 #define SEP_NOISE_VAR 2 0030 0031 /* input flags for aperture photometry */ 0032 #define SEP_MASK_IGNORE 0x0004 0033 0034 /* threshold interpretation for sep_extract */ 0035 #define SEP_THRESH_REL 0 /* in units of standard deviations (sigma) */ 0036 #define SEP_THRESH_ABS 1 /* absolute data values */ 0037 0038 /* filter types for sep_extract */ 0039 #define SEP_FILTER_CONV 0 0040 #define SEP_FILTER_MATCHED 1 0041 0042 /* structs ------------------------------------------------------------------*/ 0043 0044 /* sep_image 0045 * 0046 * Represents an image, including data, noise and mask arrays, and 0047 * gain. 0048 */ 0049 typedef struct { 0050 void *data; /* data array */ 0051 void *noise; /* noise array (can be NULL) */ 0052 void *mask; /* mask array (can be NULL) */ 0053 int dtype; /* element type of image */ 0054 int ndtype; /* element type of noise */ 0055 int mdtype; /* element type of mask */ 0056 int w; /* array width */ 0057 int h; /* array height */ 0058 double noiseval; /* scalar noise value; used only if noise == NULL */ 0059 short noise_type; /* interpretation of noise value */ 0060 double gain; /* (poisson counts / data unit) */ 0061 double maskthresh; /* pixel considered masked if mask > maskthresh */ 0062 } sep_image; 0063 0064 /* sep_bkg 0065 * 0066 * The result of sep_background() -- represents a smooth image background 0067 * and its noise with splines. 0068 */ 0069 typedef struct { 0070 int w, h; /* original image width, height */ 0071 int bw, bh; /* single tile width, height */ 0072 int nx, ny; /* number of tiles in x, y */ 0073 int n; /* nx*ny */ 0074 float global; /* global mean */ 0075 float globalrms; /* global sigma */ 0076 float *back; /* node data for interpolation */ 0077 float *dback; 0078 float *sigma; 0079 float *dsigma; 0080 } sep_bkg; 0081 0082 /* sep_catalog 0083 * 0084 * The result of sep_extract(). This is a struct of arrays. Each array has 0085 * one entry per detected object. 0086 */ 0087 typedef struct { 0088 int nobj; /* number of objects (length of all arrays) */ 0089 float *thresh; /* threshold (ADU) */ 0090 int *npix; /* # pixels extracted (size of pix array) */ 0091 int *tnpix; /* # pixels above thresh (unconvolved) */ 0092 int *xmin, *xmax; 0093 int *ymin, *ymax; 0094 double *x, *y; /* barycenter (first moments) */ 0095 double *x2, *y2, *xy; /* second moments */ 0096 double *errx2, *erry2, *errxy; /* second moment errors */ 0097 float *a, *b, *theta; /* ellipse parameters */ 0098 float *cxx, *cyy, *cxy; /* ellipse parameters (alternative) */ 0099 float *cflux; /* total flux of pixels (convolved im) */ 0100 float *flux; /* total flux of pixels (unconvolved) */ 0101 float *cpeak; /* peak intensity (ADU) (convolved) */ 0102 float *peak; /* peak intensity (ADU) (unconvolved) */ 0103 int *xcpeak, *ycpeak; /* x, y coords of peak (convolved) pixel */ 0104 int *xpeak, *ypeak; /* x, y coords of peak (unconvolved) pixel */ 0105 short *flag; /* extraction flags */ 0106 int **pix; /* array giving indicies of object's pixels in */ 0107 /* image (linearly indexed). Length is `npix`. */ 0108 /* (pointer to within the `objectspix` buffer) */ 0109 int *objectspix; /* buffer holding pixel indicies for all objects */ 0110 } sep_catalog; 0111 0112 #ifdef __cplusplus 0113 extern "C" { 0114 #endif 0115 0116 /*--------------------- global background estimation ------------------------*/ 0117 0118 /* sep_background() 0119 * 0120 * Create representation of spatially varying image background and variance. 0121 * 0122 * Note that the returned pointer must eventually be freed by calling 0123 * `sep_bkg_free()`. 0124 * 0125 * In addition to the image mask (if present), pixels <= -1e30 and NaN 0126 * are ignored. 0127 * 0128 * Source Extractor defaults: 0129 * 0130 * - bw, bh = (64, 64) 0131 * - fw, fh = (3, 3) 0132 * - fthresh = 0.0 0133 */ 0134 int sep_background(sep_image *image, 0135 int bw, int bh, /* size of a single background tile */ 0136 int fw, int fh, /* filter size in tiles */ 0137 double fthresh, /* filter threshold */ 0138 sep_bkg **bkg); /* OUTPUT */ 0139 0140 0141 /* sep_bkg_global[rms]() 0142 * 0143 * Get the estimate of the global background "median" or standard deviation. 0144 */ 0145 float sep_bkg_global(sep_bkg *bkg); 0146 float sep_bkg_globalrms(sep_bkg *bkg); 0147 0148 0149 /* sep_bkg_pix() 0150 * 0151 * Return background at (x, y). 0152 * Unlike other routines, this uses simple linear interpolation. 0153 */ 0154 float sep_bkg_pix(sep_bkg *bkg, int x, int y); 0155 0156 0157 /* sep_bkg_[sub,rms]line() 0158 * 0159 * Evaluate the background or RMS at line `y`. 0160 * Uses bicubic spline interpolation between background map vertices. 0161 * The second function subtracts the background from the input array. 0162 * Line must be an array with same width as original image. 0163 */ 0164 int sep_bkg_line(sep_bkg *bkg, int y, void *line, int dtype); 0165 int sep_bkg_subline(sep_bkg *bkg, int y, void *line, int dtype); 0166 int sep_bkg_rmsline(sep_bkg *bkg, int y, void *line, int dtype); 0167 0168 0169 /* sep_bkg_[sub,rms]array() 0170 * 0171 * Evaluate the background or RMS for entire image. 0172 * Uses bicubic spline interpolation between background map vertices. 0173 * The second function subtracts the background from the input array. 0174 * `arr` must be an array of the same size as original image. 0175 */ 0176 int sep_bkg_array(sep_bkg *bkg, void *arr, int dtype); 0177 int sep_bkg_subarray(sep_bkg *bkg, void *arr, int dtype); 0178 int sep_bkg_rmsarray(sep_bkg *bkg, void *arr, int dtype); 0179 0180 /* sep_bkg_free() 0181 * 0182 * Free memory associated with bkg. 0183 */ 0184 void sep_bkg_free(sep_bkg *bkg); 0185 0186 /*-------------------------- source extraction ------------------------------*/ 0187 0188 /* sep_extract() 0189 * 0190 * Extract sources from an image. Source Extractor defaults are shown 0191 * in [ ] above. 0192 * 0193 * Notes 0194 * ----- 0195 * `dtype` and `ndtype` indicate the data type (float, int, double) of the 0196 * image and noise arrays, respectively. 0197 * 0198 * If `noise` is NULL, thresh is interpreted as an absolute threshold. 0199 * If `noise` is not null, thresh is interpreted as a relative threshold 0200 * (the absolute threshold will be thresh*noise[i,j]). 0201 * 0202 */ 0203 int sep_extract(sep_image *image, 0204 float thresh, /* detection threshold [1.5] */ 0205 int thresh_type, /* threshold units [SEP_THRESH_REL] */ 0206 int minarea, /* minimum area in pixels [5] */ 0207 float *conv, /* convolution array (can be NULL) */ 0208 /* [{1 2 1 2 4 2 1 2 1}] */ 0209 int convw, int convh, /* w, h of convolution array [3,3] */ 0210 int filter_type, /* convolution (0) or matched (1) [0] */ 0211 int deblend_nthresh, /* deblending thresholds [32] */ 0212 double deblend_cont, /* min. deblending contrast [0.005] */ 0213 int clean_flag, /* perform cleaning? [1] */ 0214 double clean_param, /* clean parameter [1.0] */ 0215 sep_catalog **catalog); /* OUTPUT catalog */ 0216 0217 0218 0219 /* set and get the size of the pixel stack used in extract() */ 0220 void sep_set_extract_pixstack(size_t val); 0221 size_t sep_get_extract_pixstack(void); 0222 0223 /* free memory associated with a catalog */ 0224 void sep_catalog_free(sep_catalog *catalog); 0225 0226 /*-------------------------- aperture photometry ----------------------------*/ 0227 0228 0229 /* Sum array values within a circular aperture. 0230 * 0231 * Notes 0232 * ----- 0233 * error : Can be a scalar (default), an array, or NULL 0234 * If an array, set the flag SEP_ERROR_IS_ARRAY in `inflags`. 0235 * Can represent 1-sigma std. deviation (default) or variance. 0236 * If variance, set the flag SEP_ERROR_IS_VARIANCE in `inflags`. 0237 * 0238 * gain : If 0.0, poisson noise on sum is ignored when calculating error. 0239 * Otherwise, (sum / gain) is added to the variance on sum. 0240 * 0241 * area : Total pixel area included in sum. Includes masked pixels that were 0242 * corrected. The area can differ from the exact area of a circle due 0243 * to inexact subpixel sampling and intersection with array boundaries. 0244 */ 0245 int sep_sum_circle(sep_image *image, 0246 double x, /* center of aperture in x */ 0247 double y, /* center of aperture in y */ 0248 double r, /* radius of aperture */ 0249 int subpix, /* subpixel sampling */ 0250 short inflags, /* input flags (see below) */ 0251 double *sum, /* OUTPUT: sum */ 0252 double *sumerr, /* OUTPUT: error on sum */ 0253 double *area, /* OUTPUT: area included in sum */ 0254 short *flag); /* OUTPUT: flags */ 0255 0256 0257 int sep_sum_circann(sep_image *image, 0258 double x, double y, double rin, double rout, 0259 int subpix, short inflags, 0260 double *sum, double *sumerr, double *area, short *flag); 0261 0262 int sep_sum_ellipse(sep_image *image, 0263 double x, double y, double a, double b, double theta, 0264 double r, int subpix, short inflags, 0265 double *sum, double *sumerr, double *area, short *flag); 0266 0267 int sep_sum_ellipann(sep_image *image, 0268 double x, double y, double a, double b, double theta, 0269 double rin, double rout, int subpix, short inflags, 0270 double *sum, double *sumerr, double *area, short *flag); 0271 0272 /* sep_sum_circann_multi() 0273 * 0274 * Sum an array of circular annuli more efficiently (but with no exact mode). 0275 * 0276 * Notable parameters: 0277 * 0278 * rmax: Input radii are [rmax/n, 2*rmax/n, 3*rmax/n, ..., rmax]. 0279 * n: Length of input and output arrays. 0280 * sum: Preallocated array of length n holding sums in annuli. sum[0] 0281 * corresponds to r=[0, rmax/n], sum[n-1] to outermost annulus. 0282 * sumvar: Preallocated array of length n holding variance on sums. 0283 * area: Preallocated array of length n holding area summed in each annulus. 0284 * maskarea: Preallocated array of length n holding masked area in each 0285 annulus (if mask not NULL). 0286 * flag: Output flag (non-array). 0287 */ 0288 int sep_sum_circann_multi(sep_image *im, 0289 double x, double y, double rmax, int n, int subpix, 0290 short inflag, 0291 double *sum, double *sumvar, double *area, 0292 double *maskarea, short *flag); 0293 0294 /* sep_flux_radius() 0295 * 0296 * Calculate the radii enclosing the requested fraction of flux relative 0297 * to radius rmax. 0298 * 0299 * (see previous functions for most arguments) 0300 * rmax : maximum radius to analyze 0301 * fluxtot : scale requested flux fractions to this. (If NULL, flux within 0302 `rmax` is used.) 0303 * fluxfrac : array of requested fractions. 0304 * n : length of fluxfrac 0305 * r : (output) result array of length n. 0306 * flag : (output) scalar flag 0307 */ 0308 int sep_flux_radius(sep_image *im, 0309 double x, double y, double rmax, int subpix, short inflag, 0310 double *fluxtot, double *fluxfrac, int n, 0311 double *r, short *flag); 0312 0313 /* sep_kron_radius() 0314 * 0315 * Calculate Kron radius within an ellipse given by 0316 * 0317 * cxx*(x'-x)^2 + cyy*(y'-y)^2 + cxy*(x'-x)*(y'-y) < r^2 0318 * 0319 * The Kron radius is sum(r_i * v_i) / sum(v_i) where v_i is the value of pixel 0320 * i and r_i is the "radius" of pixel i, as given by the left hand side of 0321 * the above equation. 0322 * 0323 * Flags that might be set: 0324 * SEP_APER_HASMASKED - at least one of the pixels in the ellipse is masked. 0325 * SEP_APER_ALLMASKED - All pixels in the ellipse are masked. kronrad = 0. 0326 * SEP_APER_NONPOSITIVE - There was a non-positive numerator or denominator. 0327 * kronrad = 0. 0328 */ 0329 int sep_kron_radius(sep_image *im, double x, double y, 0330 double cxx, double cyy, double cxy, double r, 0331 double *kronrad, short *flag); 0332 0333 0334 /* sep_windowed() 0335 * 0336 * Calculate "windowed" position parameters via an iterative procedure. 0337 * 0338 * x, y : initial center 0339 * sig : sigma of Gaussian to use for weighting. The integration 0340 * radius is 4 * sig. 0341 * subpix : Subpixels to use in aperture-pixel overlap. 0342 * SExtractor uses 11. 0 is supported for exact overlap. 0343 * xout, yout : output center. 0344 * niter : number of iterations used. 0345 */ 0346 int sep_windowed(sep_image *im, 0347 double x, double y, double sig, int subpix, short inflag, 0348 double *xout, double *yout, int *niter, short *flag); 0349 0350 0351 /* sep_set_ellipse() 0352 * 0353 * Set array elements within an elliptical aperture to a given value. 0354 * 0355 * Ellipse: cxx*(x'-x)^2 + cyy*(y'-y)^2 + cxy*(x'-x)*(y'-y) = r^2 0356 */ 0357 void sep_set_ellipse(unsigned char *arr, int w, int h, 0358 double x, double y, double cxx, double cyy, double cxy, 0359 double r, unsigned char val); 0360 0361 0362 /* sep_ellipse_axes() 0363 * sep_ellipse_coeffs() 0364 * 0365 * Convert between coefficient representation of ellipse, 0366 * cxx*(x'-x)^2 + cyy*(y'-y)^2 + cxy*(x'-x)*(y'-y) = r^2, 0367 * and axis representation of an ellipse. The axis representation is 0368 * defined by: 0369 * 0370 * a = semimajor axis 0371 * b = semiminor axis 0372 * theta = angle in radians counter-clockwise from positive x axis 0373 */ 0374 int sep_ellipse_axes(double cxx, double cyy, double cxy, 0375 double *a, double *b, double *theta); 0376 void sep_ellipse_coeffs(double a, double b, double theta, 0377 double *cxx, double *cyy, double *cxy); 0378 0379 /*----------------------- info & error messaging ----------------------------*/ 0380 0381 /* sep_version_string : library version (e.g., "0.2.0") */ 0382 extern char *sep_version_string; 0383 0384 /* sep_get_errmsg() 0385 * 0386 * Return a short descriptive error message that corresponds to the input 0387 * error status value. The message may be up to 60 characters long, plus 0388 * the terminating null character. 0389 */ 0390 void sep_get_errmsg(int status, char *errtext); 0391 0392 0393 /* sep_get_errdetail() 0394 * 0395 * Return a longer error message with more specifics about the problem. 0396 * The message may be up to 512 characters. 0397 */ 0398 void sep_get_errdetail(char *errtext); 0399 0400 #ifdef __cplusplus 0401 } 0402 #endif