File indexing completed on 2024-04-28 03:43:53
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 #include <math.h> 0011 #include <stdio.h> 0012 #include <stdlib.h> 0013 #include <string.h> 0014 #include "sep.h" 0015 #include "sepcore.h" 0016 #include "overlap.h" 0017 0018 #define FLUX_RADIUS_BUFSIZE 64 0019 0020 #define WINPOS_NITERMAX 16 /* Maximum number of steps */ 0021 #define WINPOS_NSIG 4 /* Measurement radius */ 0022 #define WINPOS_STEPMIN 0.0001 /* Minimum change in position for continuing */ 0023 #define WINPOS_FAC 2.0 /* Centroid offset factor (2 for a Gaussian) */ 0024 0025 /* 0026 Adding (void *) pointers is a GNU C extension, not part of standard C. 0027 When compiling on Windows with MS VIsual C compiler need to cast the 0028 (void *) to something the size of one byte. 0029 */ 0030 #if defined(_MSC_VER) 0031 #define MSVC_VOID_CAST (char *) 0032 #else 0033 #define MSVC_VOID_CAST 0034 #endif 0035 0036 /****************************************************************************/ 0037 /* conversions between ellipse representations */ 0038 0039 /* return ellipse semi-major and semi-minor axes and position angle, given 0040 representation: cxx*x^2 + cyy*y^2 + cxy*x*y = 1 0041 derived from http://mathworld.wolfram.com/Ellipse.html 0042 0043 Input requirements: 0044 cxx*cyy - cxy*cxy/4. > 0. 0045 cxx + cyy > 0. 0046 */ 0047 int sep_ellipse_axes(double cxx, double cyy, double cxy, 0048 double *a, double *b, double *theta) 0049 { 0050 double p, q, t; 0051 0052 p = cxx + cyy; 0053 q = cxx - cyy; 0054 t = sqrt(q*q + cxy*cxy); 0055 0056 /* Ensure that parameters actually describe an ellipse. */ 0057 if ((cxx*cyy - cxy*cxy/4. <= 0.) || (p <= 0.)) 0058 return NON_ELLIPSE_PARAMS; 0059 0060 *a = sqrt(2. / (p - t)); 0061 *b = sqrt(2. / (p + t)); 0062 0063 /* theta = 0 if cxy == 0, else (1/2) acot(q/cxy) */ 0064 *theta = (cxy == 0.) ? 0. : (q == 0. ? 0. : atan(cxy/q))/2.; 0065 if (cxx>cyy) 0066 *theta += PI/2.; 0067 if (*theta > PI/2.) 0068 *theta -= PI; 0069 0070 return RETURN_OK; 0071 } 0072 0073 void sep_ellipse_coeffs(double a, double b, double theta, 0074 double *cxx, double *cyy, double *cxy) 0075 { 0076 double costheta, sintheta; 0077 0078 costheta = cos(theta); 0079 sintheta = sin(theta); 0080 0081 *cxx = costheta*costheta/(a*a) + sintheta*sintheta/(b*b); 0082 *cyy = sintheta*sintheta/(a*a) + costheta*costheta/(b*b); 0083 *cxy = 2.*costheta*sintheta * (1./(a*a) - 1./(b*b)); 0084 } 0085 0086 /*****************************************************************************/ 0087 /* Helper functions for aperture functions */ 0088 0089 /* determine the extent of the box that just contains the circle with 0090 * parameters x, y, r. xmin is inclusive and xmax is exclusive. 0091 * Ensures that box is within image bound and sets a flag if it is not. 0092 */ 0093 static void boxextent(double x, double y, double rx, double ry, int w, int h, 0094 int *xmin, int *xmax, int *ymin, int *ymax, 0095 short *flag) 0096 { 0097 *xmin = (int)(x - rx + 0.5); 0098 *xmax = (int)(x + rx + 1.4999999); 0099 *ymin = (int)(y - ry + 0.5); 0100 *ymax = (int)(y + ry + 1.4999999); 0101 if (*xmin < 0) 0102 { 0103 *xmin = 0; 0104 *flag |= SEP_APER_TRUNC; 0105 } 0106 if (*xmax > w) 0107 { 0108 *xmax = w; 0109 *flag |= SEP_APER_TRUNC; 0110 } 0111 if (*ymin < 0) 0112 { 0113 *ymin = 0; 0114 *flag |= SEP_APER_TRUNC; 0115 } 0116 if (*ymax > h) 0117 { 0118 *ymax = h; 0119 *flag |= SEP_APER_TRUNC; 0120 } 0121 } 0122 0123 0124 static void boxextent_ellipse(double x, double y, 0125 double cxx, double cyy, double cxy, double r, 0126 int w, int h, 0127 int *xmin, int *xmax, int *ymin, int *ymax, 0128 short *flag) 0129 { 0130 double dxlim, dylim; 0131 0132 dxlim = cxx - cxy*cxy/(4.0*cyy); 0133 dxlim = dxlim>0.0 ? r/sqrt(dxlim) : 0.0; 0134 dylim = cyy - cxy*cxy/(4.0*cxx); 0135 dylim = dylim > 0.0 ? r/sqrt(dylim) : 0.0; 0136 boxextent(x, y, dxlim, dylim, w, h, xmin, xmax, ymin, ymax, flag); 0137 } 0138 0139 /* determine oversampled annulus for a circle */ 0140 static void oversamp_ann_circle(double r, double *r_in2, double *r_out2) 0141 { 0142 *r_in2 = r - 0.7072; 0143 *r_in2 = (*r_in2 > 0.0) ? (*r_in2)*(*r_in2) : 0.0; 0144 *r_out2 = r + 0.7072; 0145 *r_out2 = (*r_out2) * (*r_out2); 0146 } 0147 0148 /* determine oversampled "annulus" for an ellipse */ 0149 static void oversamp_ann_ellipse(double r, double b, double *r_in2, 0150 double *r_out2) 0151 { 0152 *r_in2 = r - 0.7072/b; 0153 *r_in2 = (*r_in2 > 0.0) ? (*r_in2)*(*r_in2) : 0.0; 0154 *r_out2 = r + 0.7072/b; 0155 *r_out2 = (*r_out2) * (*r_out2); 0156 } 0157 0158 /*****************************************************************************/ 0159 /* circular aperture */ 0160 0161 #define APER_NAME sep_sum_circle 0162 #define APER_ARGS double r 0163 #define APER_DECL double r2, r_in2, r_out2 0164 #define APER_CHECKS \ 0165 if (r < 0.0) \ 0166 return ILLEGAL_APER_PARAMS 0167 #define APER_INIT \ 0168 r2 = r*r; \ 0169 oversamp_ann_circle(r, &r_in2, &r_out2) 0170 #define APER_BOXEXTENT boxextent(x, y, r, r, im->w, im->h, \ 0171 &xmin, &xmax, &ymin, &ymax, flag) 0172 #define APER_EXACT circoverlap(dx-0.5, dy-0.5, dx+0.5, dy+0.5, r) 0173 #define APER_RPIX2 dx*dx + dy*dy 0174 #define APER_RPIX2_SUBPIX dx1*dx1 + dy2 0175 #define APER_COMPARE1 rpix2 < r_out2 0176 #define APER_COMPARE2 rpix2 > r_in2 0177 #define APER_COMPARE3 rpix2 < r2 0178 #include "aperture.i" 0179 #undef APER_NAME 0180 #undef APER_ARGS 0181 #undef APER_DECL 0182 #undef APER_CHECKS 0183 #undef APER_INIT 0184 #undef APER_BOXEXTENT 0185 #undef APER_EXACT 0186 #undef APER_RPIX2 0187 #undef APER_RPIX2_SUBPIX 0188 #undef APER_COMPARE1 0189 #undef APER_COMPARE2 0190 #undef APER_COMPARE3 0191 0192 /*****************************************************************************/ 0193 /* elliptical aperture */ 0194 0195 #define APER_NAME sep_sum_ellipse 0196 #define APER_ARGS double a, double b, double theta, double r 0197 #define APER_DECL double cxx, cyy, cxy, r2, r_in2, r_out2 0198 #define APER_CHECKS \ 0199 if (!(r >= 0.0 && b >= 0.0 && a >= b && \ 0200 theta >= -PI/2. && theta <= PI/2.)) \ 0201 return ILLEGAL_APER_PARAMS 0202 #define APER_INIT \ 0203 r2 = r*r; \ 0204 oversamp_ann_ellipse(r, b, &r_in2, &r_out2); \ 0205 sep_ellipse_coeffs(a, b, theta, &cxx, &cyy, &cxy); \ 0206 a *= r; \ 0207 b *= r 0208 #define APER_BOXEXTENT boxextent_ellipse(x, y, cxx, cyy, cxy, r, im->w, im->h, \ 0209 &xmin, &xmax, &ymin, &ymax, flag) 0210 #define APER_EXACT ellipoverlap(dx-0.5, dy-0.5, dx+0.5, dy+0.5, a, b, theta) 0211 #define APER_RPIX2 cxx*dx*dx + cyy*dy*dy + cxy*dx*dy 0212 #define APER_RPIX2_SUBPIX cxx*dx1*dx1 + cyy*dy2 + cxy*dx1*dy 0213 #define APER_COMPARE1 rpix2 < r_out2 0214 #define APER_COMPARE2 rpix2 > r_in2 0215 #define APER_COMPARE3 rpix2 < r2 0216 #include "aperture.i" 0217 #undef APER_NAME 0218 #undef APER_ARGS 0219 #undef APER_DECL 0220 #undef APER_CHECKS 0221 #undef APER_INIT 0222 #undef APER_BOXEXTENT 0223 #undef APER_EXACT 0224 #undef APER_RPIX2 0225 #undef APER_RPIX2_SUBPIX 0226 #undef APER_COMPARE1 0227 #undef APER_COMPARE2 0228 #undef APER_COMPARE3 0229 0230 /*****************************************************************************/ 0231 /* circular annulus aperture */ 0232 0233 #define APER_NAME sep_sum_circann 0234 #define APER_ARGS double rin, double rout 0235 #define APER_DECL double rin2, rin_in2, rin_out2, rout2, rout_in2, rout_out2 0236 #define APER_CHECKS \ 0237 if (!(rin >= 0.0 && rout >= rin)) \ 0238 return ILLEGAL_APER_PARAMS 0239 #define APER_INIT \ 0240 rin2 = rin*rin; \ 0241 oversamp_ann_circle(rin, &rin_in2, &rin_out2); \ 0242 rout2 = rout*rout; \ 0243 oversamp_ann_circle(rout, &rout_in2, &rout_out2) 0244 #define APER_BOXEXTENT boxextent(x, y, rout, rout, im->w, im->h, \ 0245 &xmin, &xmax, &ymin, &ymax, flag) 0246 #define APER_EXACT (circoverlap(dx-0.5, dy-0.5, dx+0.5, dy+0.5, rout) - \ 0247 circoverlap(dx-0.5, dy-0.5, dx+0.5, dy+0.5, rin)) 0248 #define APER_RPIX2 dx*dx + dy*dy 0249 #define APER_RPIX2_SUBPIX dx1*dx1 + dy2 0250 #define APER_COMPARE1 (rpix2 < rout_out2) && (rpix2 > rin_in2) 0251 #define APER_COMPARE2 (rpix2 > rout_in2) || (rpix2 < rin_out2) 0252 #define APER_COMPARE3 (rpix2 < rout2) && (rpix2 > rin2) 0253 #include "aperture.i" 0254 #undef APER_NAME 0255 #undef APER_ARGS 0256 #undef APER_DECL 0257 #undef APER_CHECKS 0258 #undef APER_INIT 0259 #undef APER_BOXEXTENT 0260 #undef APER_EXACT 0261 #undef APER_RPIX2 0262 #undef APER_RPIX2_SUBPIX 0263 #undef APER_COMPARE1 0264 #undef APER_COMPARE2 0265 #undef APER_COMPARE3 0266 0267 /*****************************************************************************/ 0268 /* elliptical annulus aperture */ 0269 0270 #define APER_NAME sep_sum_ellipann 0271 #define APER_ARGS double a, double b, double theta, double rin, double rout 0272 #define APER_DECL \ 0273 double cxx, cyy, cxy; \ 0274 double rin2, rin_in2, rin_out2, rout2, rout_in2, rout_out2 0275 #define APER_CHECKS \ 0276 if (!(rin >= 0.0 && rout >= rin && b >= 0.0 && a >= b && \ 0277 theta >= -PI/2. && theta <= PI/2.)) \ 0278 return ILLEGAL_APER_PARAMS 0279 #define APER_INIT \ 0280 rin2 = rin*rin; \ 0281 oversamp_ann_ellipse(rin, b, &rin_in2, &rin_out2); \ 0282 rout2 = rout*rout; \ 0283 oversamp_ann_ellipse(rout, b, &rout_in2, &rout_out2); \ 0284 sep_ellipse_coeffs(a, b, theta, &cxx, &cyy, &cxy) 0285 #define APER_BOXEXTENT boxextent_ellipse(x, y, cxx, cyy, cxy, rout, im->w, im->h, \ 0286 &xmin, &xmax, &ymin, &ymax, flag) 0287 #define APER_EXACT \ 0288 (ellipoverlap(dx-0.5, dy-0.5, dx+0.5, dy+0.5, a*rout, b*rout, theta) - \ 0289 ellipoverlap(dx-0.5, dy-0.5, dx+0.5, dy+0.5, a*rin, b*rin, theta)) 0290 #define APER_RPIX2 cxx*dx*dx + cyy*dy*dy + cxy*dx*dy 0291 #define APER_RPIX2_SUBPIX cxx*dx1*dx1 + cyy*dy2 + cxy*dx1*dy 0292 #define APER_COMPARE1 (rpix2 < rout_out2) && (rpix2 > rin_in2) 0293 #define APER_COMPARE2 (rpix2 > rout_in2) || (rpix2 < rin_out2) 0294 #define APER_COMPARE3 (rpix2 < rout2) && (rpix2 > rin2) 0295 #include "aperture.i" 0296 #undef APER_NAME 0297 #undef APER_ARGS 0298 #undef APER_DECL 0299 #undef APER_CHECKS 0300 #undef APER_INIT 0301 #undef APER_BOXEXTENT 0302 #undef APER_EXACT 0303 #undef APER_RPIX2 0304 #undef APER_RPIX2_SUBPIX 0305 #undef APER_COMPARE1 0306 #undef APER_COMPARE2 0307 #undef APER_COMPARE3 0308 0309 0310 /*****************************************************************************/ 0311 /* 0312 * This is just different enough from the other aperture functions 0313 * that it doesn't quite make sense to use aperture.i. 0314 */ 0315 int sep_sum_circann_multi(sep_image *im, 0316 double x, double y, double rmax, int n, int subpix, 0317 short inflag, 0318 double *sum, double *sumvar, double *area, 0319 double *maskarea, short *flag) 0320 { 0321 PIXTYPE pix, varpix; 0322 double dx, dy, dx1, dy2, offset, scale, scale2, tmp, rpix2; 0323 int ix, iy, xmin, xmax, ymin, ymax, sx, sy, status, size, esize, msize; 0324 long pos; 0325 short errisarray, errisstd; 0326 BYTE *datat, *errort, *maskt; 0327 converter convert, econvert, mconvert; 0328 double rpix, r_out, r_out2, d, prevbinmargin, nextbinmargin, step, stepdens; 0329 int j, ismasked; 0330 0331 /* input checks */ 0332 if (rmax < 0.0 || n < 1) 0333 return ILLEGAL_APER_PARAMS; 0334 if (subpix < 1) 0335 return ILLEGAL_SUBPIX; 0336 0337 /* clear results arrays */ 0338 memset(sum, 0, (size_t)(n*sizeof(double))); 0339 memset(sumvar, 0, (size_t)(n*sizeof(double))); 0340 memset(area, 0, (size_t)(n*sizeof(double))); 0341 if (im->mask) 0342 memset(maskarea, 0, (size_t)(n*sizeof(double))); 0343 0344 /* initializations */ 0345 size = esize = msize = 0; 0346 datat = maskt = NULL; 0347 errort = im->noise; 0348 *flag = 0; 0349 varpix = 0.0; 0350 scale = 1.0/subpix; 0351 scale2 = scale*scale; 0352 offset = 0.5*(scale-1.0); 0353 0354 r_out = rmax + 1.5; /* margin for interpolation */ 0355 r_out2 = r_out * r_out; 0356 step = rmax/n; 0357 stepdens = 1.0/step; 0358 prevbinmargin = 0.7072; 0359 nextbinmargin = step - 0.7072; 0360 j = 0; 0361 d = 0.; 0362 ismasked = 0; 0363 errisarray = 0; 0364 errisstd = 0; 0365 0366 /* get data converter(s) for input array(s) */ 0367 if ((status = get_converter(im->dtype, &convert, &size))) 0368 return status; 0369 if (im->mask && (status = get_converter(im->mdtype, &mconvert, &msize))) 0370 return status; 0371 0372 /* get image noise */ 0373 if (im->noise_type != SEP_NOISE_NONE) 0374 { 0375 errisstd = (im->noise_type == SEP_NOISE_STDDEV); 0376 if (im->noise) 0377 { 0378 errisarray = 1; 0379 if ((status = get_converter(im->ndtype, &econvert, &esize))) 0380 return status; 0381 } 0382 else 0383 { 0384 varpix = (errisstd)? im->noiseval * im->noiseval: im->noiseval; 0385 } 0386 } 0387 0388 0389 /* get extent of box */ 0390 boxextent(x, y, r_out, r_out, im->w, im->h, &xmin, &xmax, &ymin, &ymax, 0391 flag); 0392 0393 /* loop over rows in the box */ 0394 for (iy=ymin; iy<ymax; iy++) 0395 { 0396 /* set pointers to the start of this row */ 0397 pos = (iy % im->h) * im->w + xmin; 0398 datat = MSVC_VOID_CAST im->data + pos*size; 0399 if (errisarray) 0400 errort = MSVC_VOID_CAST im->noise + pos*esize; 0401 if (im->mask) 0402 maskt = MSVC_VOID_CAST im->mask + pos*msize; 0403 0404 /* loop over pixels in this row */ 0405 for (ix=xmin; ix<xmax; ix++) 0406 { 0407 dx = ix - x; 0408 dy = iy - y; 0409 rpix2 = dx*dx + dy*dy; 0410 if (rpix2 < r_out2) 0411 { 0412 /* get pixel values */ 0413 pix = convert(datat); 0414 if (errisarray) 0415 { 0416 varpix = econvert(errort); 0417 if (errisstd) 0418 varpix *= varpix; 0419 } 0420 if (im->mask) 0421 { 0422 if (mconvert(maskt) > im->maskthresh) 0423 { 0424 *flag |= SEP_APER_HASMASKED; 0425 ismasked = 1; 0426 } 0427 else 0428 ismasked = 0; 0429 } 0430 0431 /* check if oversampling is needed (close to bin boundary?) */ 0432 rpix = sqrt(rpix2); 0433 d = fmod(rpix, step); 0434 if (d < prevbinmargin || d > nextbinmargin) 0435 { 0436 dx += offset; 0437 dy += offset; 0438 for (sy=subpix; sy--; dy+=scale) 0439 { 0440 dx1 = dx; 0441 dy2 = dy*dy; 0442 for (sx=subpix; sx--; dx1+=scale) 0443 { 0444 j = (int)(sqrt(dx1*dx1+dy2)*stepdens); 0445 if (j < n) 0446 { 0447 if (ismasked) 0448 maskarea[j] += scale2; 0449 else 0450 { 0451 sum[j] += scale2*pix; 0452 sumvar[j] += scale2*varpix; 0453 } 0454 area[j] += scale2; 0455 } 0456 } 0457 } 0458 } 0459 else 0460 /* pixel not close to bin boundary */ 0461 { 0462 j = (int)(rpix*stepdens); 0463 if (j < n) 0464 { 0465 if (ismasked) 0466 maskarea[j] += 1.0; 0467 else 0468 { 0469 sum[j] += pix; 0470 sumvar[j] += varpix; 0471 } 0472 area[j] += 1.0; 0473 } 0474 } 0475 } /* closes "if pixel might be within aperture" */ 0476 0477 /* increment pointers by one element */ 0478 datat += size; 0479 if (errisarray) 0480 errort += esize; 0481 maskt += msize; 0482 } 0483 } 0484 0485 0486 /* correct for masked values */ 0487 if (im->mask) 0488 { 0489 if (inflag & SEP_MASK_IGNORE) 0490 for (j=n; j--;) 0491 area[j] -= maskarea[j]; 0492 else 0493 { 0494 for (j=n; j--;) 0495 { 0496 tmp = area[j] == maskarea[j]? 0.0: area[j]/(area[j]-maskarea[j]); 0497 sum[j] *= tmp; 0498 sumvar[j] *= tmp; 0499 } 0500 } 0501 } 0502 0503 /* add poisson noise, only if gain > 0 */ 0504 if (im->gain > 0.0) 0505 for (j=n; j--;) 0506 if (sum[j] > 0.0) 0507 sumvar[j] += sum[j] / im->gain; 0508 0509 return status; 0510 } 0511 0512 0513 /* for use in flux_radius */ 0514 static double inverse(double xmax, double *y, int n, double ytarg) 0515 { 0516 double step; 0517 int i; 0518 0519 step = xmax/n; 0520 i = 0; 0521 0522 /* increment i until y[i] is >= to ytarg */ 0523 while (i < n && y[i] < ytarg) 0524 i++; 0525 0526 if (i == 0) 0527 { 0528 if (ytarg <= 0. || y[0] == 0.) 0529 return 0.; 0530 return step * ytarg/y[0]; 0531 } 0532 if (i == n) 0533 return xmax; 0534 0535 /* note that y[i-1] corresponds to x=step*i. */ 0536 return step * (i + (ytarg - y[i-1])/(y[i] - y[i-1])); 0537 } 0538 0539 int sep_flux_radius(sep_image *im, 0540 double x, double y, double rmax, int subpix, short inflag, 0541 double *fluxtot, double *fluxfrac, int n, double *r, 0542 short *flag) 0543 { 0544 int status; 0545 int i; 0546 double f; 0547 double sumbuf[FLUX_RADIUS_BUFSIZE] = {0.}; 0548 double sumvarbuf[FLUX_RADIUS_BUFSIZE]; 0549 double areabuf[FLUX_RADIUS_BUFSIZE]; 0550 double maskareabuf[FLUX_RADIUS_BUFSIZE]; 0551 0552 /* measure FLUX_RADIUS_BUFSIZE annuli out to rmax. */ 0553 status = sep_sum_circann_multi(im, x, y, rmax, 0554 FLUX_RADIUS_BUFSIZE, subpix, inflag, 0555 sumbuf, sumvarbuf, areabuf, maskareabuf, 0556 flag); 0557 0558 /* sum up sumbuf array */ 0559 for (i=1; i<FLUX_RADIUS_BUFSIZE; i++) 0560 sumbuf[i] += sumbuf[i-1]; 0561 0562 /* if given, use "total flux", else, use sum within rmax. */ 0563 f = fluxtot? *fluxtot : sumbuf[FLUX_RADIUS_BUFSIZE-1]; 0564 0565 /* Use inverse to get the radii corresponding to the requested flux fracs */ 0566 for (i=0; i<n; i++) 0567 r[i] = inverse(rmax, sumbuf, FLUX_RADIUS_BUFSIZE, fluxfrac[i] * f); 0568 0569 return status; 0570 } 0571 0572 /*****************************************************************************/ 0573 /* calculate Kron radius from pixels within an ellipse. */ 0574 int sep_kron_radius(sep_image *im, double x, double y, 0575 double cxx, double cyy, double cxy, double r, 0576 double *kronrad, short *flag) 0577 { 0578 float pix; 0579 double r1, v1, r2, area, rpix2, dx, dy; 0580 int ix, iy, xmin, xmax, ymin, ymax, status, size, msize; 0581 long pos; 0582 BYTE *datat, *maskt; 0583 converter convert, mconvert; 0584 0585 r2 = r*r; 0586 r1 = v1 = 0.0; 0587 area = 0.0; 0588 *flag = 0; 0589 datat = maskt = NULL; 0590 size = msize = 0; 0591 0592 /* get data converter(s) for input array(s) */ 0593 if ((status = get_converter(im->dtype, &convert, &size))) 0594 return status; 0595 if (im->mask && (status = get_converter(im->mdtype, &mconvert, &msize))) 0596 return status; 0597 0598 0599 /* get extent of ellipse in x and y */ 0600 boxextent_ellipse(x, y, cxx, cyy, cxy, r, im->w, im->h, 0601 &xmin, &xmax, &ymin, &ymax, flag); 0602 0603 /* loop over rows in the box */ 0604 for (iy=ymin; iy<ymax; iy++) 0605 { 0606 /* set pointers to the start of this row */ 0607 pos = (iy % im->h) * im->w + xmin; 0608 datat = MSVC_VOID_CAST im->data + pos*size; 0609 if (im->mask) 0610 maskt = MSVC_VOID_CAST im->mask + pos*msize; 0611 0612 /* loop over pixels in this row */ 0613 for (ix=xmin; ix<xmax; ix++) 0614 { 0615 dx = ix - x; 0616 dy = iy - y; 0617 rpix2 = cxx*dx*dx + cyy*dy*dy + cxy*dx*dy; 0618 if (rpix2 <= r2) 0619 { 0620 pix = convert(datat); 0621 if ((pix < -BIG) || (im->mask && mconvert(maskt) > im->maskthresh)) 0622 { 0623 *flag |= SEP_APER_HASMASKED; 0624 } 0625 else 0626 { 0627 r1 += sqrt(rpix2)*pix; 0628 v1 += pix; 0629 area++; 0630 } 0631 } 0632 0633 /* increment pointers by one element */ 0634 datat += size; 0635 maskt += msize; 0636 } 0637 } 0638 0639 if (area == 0) 0640 { 0641 *flag |= SEP_APER_ALLMASKED; 0642 *kronrad = 0.0; 0643 } 0644 else if (r1 <= 0.0 || v1 <= 0.0) 0645 { 0646 *flag |= SEP_APER_NONPOSITIVE; 0647 *kronrad = 0.0; 0648 } 0649 else 0650 { 0651 *kronrad = r1 / v1; 0652 } 0653 0654 return RETURN_OK; 0655 } 0656 0657 0658 /* set array values within an ellipse (uc = unsigned char array) */ 0659 void sep_set_ellipse(unsigned char *arr, int w, int h, 0660 double x, double y, double cxx, double cyy, double cxy, 0661 double r, unsigned char val) 0662 { 0663 unsigned char *arrt; 0664 int xmin, xmax, ymin, ymax, xi, yi; 0665 double r2, dx, dy, dy2; 0666 short flag; /* not actually used, but input to boxextent */ 0667 0668 flag = 0; 0669 r2 = r*r; 0670 0671 boxextent_ellipse(x, y, cxx, cyy, cxy, r, w, h, 0672 &xmin, &xmax, &ymin, &ymax, &flag); 0673 0674 for (yi=ymin; yi<ymax; yi++) 0675 { 0676 arrt = arr + (yi*w + xmin); 0677 dy = yi - y; 0678 dy2 = dy*dy; 0679 for (xi=xmin; xi<xmax; xi++, arrt++) 0680 { 0681 dx = xi - x; 0682 if ((cxx*dx*dx + cyy*dy2 + cxy*dx*dy) <= r2) 0683 *arrt = val; 0684 } 0685 } 0686 } 0687 0688 /*****************************************************************************/ 0689 /* 0690 * As with `sep_sum_circann_multi`, this is different enough from the other 0691 * aperture functions that it doesn't quite make sense to use aperture.i. 0692 * 0693 * To reproduce SExtractor: 0694 * - use sig = obj.hl_radius * 2/2.35. 0695 * - use obj.posx/posy for initial position 0696 * 0697 */ 0698 0699 int sep_windowed(sep_image *im, 0700 double x, double y, double sig, int subpix, short inflag, 0701 double *xout, double *yout, int *niter, short *flag) 0702 { 0703 PIXTYPE pix, varpix; 0704 double dx, dy, dx1, dy2, offset, scale, scale2, tmp, dxpos, dypos, weight; 0705 double maskarea, maskweight, maskdxpos, maskdypos; 0706 double r, tv, twv, sigtv, totarea, overlap, rpix2, invtwosig2; 0707 double wpix; 0708 int i, ix, iy, xmin, xmax, ymin, ymax, sx, sy, status, size, esize, msize; 0709 long pos; 0710 short errisarray, errisstd; 0711 BYTE *datat, *errort, *maskt; 0712 converter convert, econvert, mconvert; 0713 double r2, r_in2, r_out2; 0714 0715 /* input checks */ 0716 if (sig < 0.0) 0717 return ILLEGAL_APER_PARAMS; 0718 if (subpix < 0) 0719 return ILLEGAL_SUBPIX; 0720 0721 /* initializations */ 0722 size = esize = msize = 0; 0723 tv = sigtv = 0.0; 0724 overlap = totarea = maskweight = 0.0; 0725 datat = maskt = NULL; 0726 errort = im->noise; 0727 *flag = 0; 0728 varpix = 0.0; 0729 scale = 1.0/subpix; 0730 scale2 = scale*scale; 0731 offset = 0.5*(scale-1.0); 0732 invtwosig2 = 1.0/(2.0*sig*sig); 0733 errisarray = 0; 0734 errisstd = 0; 0735 0736 /* Integration radius */ 0737 r = WINPOS_NSIG*sig; 0738 0739 /* calculate oversampled annulus */ 0740 r2 = r*r; 0741 oversamp_ann_circle(r, &r_in2, &r_out2); 0742 0743 /* get data converter(s) for input array(s) */ 0744 if ((status = get_converter(im->dtype, &convert, &size))) 0745 return status; 0746 if (im->mask && (status = get_converter(im->mdtype, &mconvert, &msize))) 0747 return status; 0748 0749 /* get image noise */ 0750 if (im->noise_type != SEP_NOISE_NONE) 0751 { 0752 errisstd = (im->noise_type == SEP_NOISE_STDDEV); 0753 if (im->noise) 0754 { 0755 errisarray = 1; 0756 if ((status = get_converter(im->ndtype, &econvert, &esize))) 0757 return status; 0758 } 0759 else 0760 { 0761 varpix = (errisstd)? im->noiseval * im->noiseval: im->noiseval; 0762 } 0763 } 0764 0765 /* iteration loop */ 0766 for (i=0; i<WINPOS_NITERMAX; i++) 0767 { 0768 0769 /* get extent of box */ 0770 boxextent(x, y, r, r, im->w, im->h, 0771 &xmin, &xmax, &ymin, &ymax, flag); 0772 0773 /* TODO: initialize values */ 0774 /* mx2ph 0775 my2ph 0776 esum, emxy, emx2, emy2, mx2, my2, mxy */ 0777 tv = twv = sigtv = 0.0; 0778 overlap = totarea = maskarea = maskweight = 0.0; 0779 dxpos = dypos = 0.0; 0780 maskdxpos = maskdypos = 0.0; 0781 0782 /* loop over rows in the box */ 0783 for (iy=ymin; iy<ymax; iy++) 0784 { 0785 /* set pointers to the start of this row */ 0786 pos = (iy % im->h) * im->w + xmin; 0787 datat = MSVC_VOID_CAST im->data + pos*size; 0788 if (errisarray) 0789 errort = MSVC_VOID_CAST im->noise + pos*esize; 0790 if (im->mask) 0791 maskt = MSVC_VOID_CAST im->mask + pos*msize; 0792 0793 /* loop over pixels in this row */ 0794 for (ix=xmin; ix<xmax; ix++) 0795 { 0796 dx = ix - x; 0797 dy = iy - y; 0798 rpix2 = dx*dx + dy*dy; 0799 if (rpix2 < r_out2) 0800 { 0801 if (rpix2 > r_in2) /* might be partially in aperture */ 0802 { 0803 if (subpix == 0) 0804 overlap = circoverlap(dx-0.5, dy-0.5, dx+0.5, dy+0.5, 0805 r); 0806 else 0807 { 0808 dx += offset; 0809 dy += offset; 0810 overlap = 0.0; 0811 for (sy=subpix; sy--; dy+=scale) 0812 { 0813 dx1 = dx; 0814 dy2 = dy*dy; 0815 for (sx=subpix; sx--; dx1+=scale) 0816 if (dx1*dx1 + dy2 < r2) 0817 overlap += scale2; 0818 } 0819 } 0820 } 0821 else 0822 /* definitely fully in aperture */ 0823 overlap = 1.0; 0824 0825 /* get pixel value and variance value */ 0826 pix = convert(datat); 0827 if (errisarray) 0828 { 0829 varpix = econvert(errort); 0830 if (errisstd) 0831 varpix *= varpix; 0832 } 0833 0834 /* offset of this pixel from center */ 0835 dx = ix - x; 0836 dy = iy - y; 0837 0838 /* weight by gaussian */ 0839 weight = exp(-rpix2*invtwosig2); 0840 0841 if (im->mask && (mconvert(maskt) > im->maskthresh)) 0842 { 0843 *flag |= SEP_APER_HASMASKED; 0844 maskarea += overlap; 0845 maskweight += overlap * weight; 0846 maskdxpos += overlap * weight * dx; 0847 maskdypos += overlap * weight * dy; 0848 } 0849 else 0850 { 0851 tv += pix * overlap; 0852 wpix = pix * overlap * weight; 0853 twv += wpix; 0854 dxpos += wpix * dx; 0855 dypos += wpix * dy; 0856 } 0857 0858 totarea += overlap; 0859 0860 } /* closes "if pixel might be within aperture" */ 0861 0862 /* increment pointers by one element */ 0863 datat += size; 0864 if (errisarray) 0865 errort += esize; 0866 maskt += msize; 0867 } /* closes loop over x */ 0868 } /* closes loop over y */ 0869 0870 /* we're done looping over pixels for this iteration. 0871 * Our summary statistics are: 0872 * 0873 * tv : total value 0874 * twv : total weighted value 0875 * dxpos : weighted dx 0876 * dypos : weighted dy 0877 */ 0878 0879 /* Correct for masked values: This effectively makes it as if 0880 * the masked pixels had the value of the average unmasked value 0881 * in the aperture. 0882 */ 0883 if (im->mask) 0884 { 0885 /* this option will probably not yield accurate values */ 0886 if (inflag & SEP_MASK_IGNORE) 0887 totarea -= maskarea; 0888 else 0889 { 0890 tmp = tv/(totarea-maskarea); /* avg unmasked pixel value */ 0891 twv += tmp * maskweight; 0892 dxpos += tmp * maskdxpos; 0893 dypos += tmp * maskdypos; 0894 } 0895 } 0896 0897 /* update center */ 0898 if (twv>0.0) 0899 { 0900 x += (dxpos /= twv) * WINPOS_FAC; 0901 y += (dypos /= twv) * WINPOS_FAC; 0902 } 0903 else 0904 break; 0905 0906 /* Stop here if position does not change */ 0907 if (dxpos*dxpos+dypos*dypos < WINPOS_STEPMIN*WINPOS_STEPMIN) 0908 break; 0909 0910 } /* closes loop over interations */ 0911 0912 /* assign output results */ 0913 *xout = x; 0914 *yout = y; 0915 *niter = i+1; 0916 0917 return status; 0918 }