From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- noao/nproto/ace/evaluate.x | 641 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 641 insertions(+) create mode 100644 noao/nproto/ace/evaluate.x (limited to 'noao/nproto/ace/evaluate.x') diff --git a/noao/nproto/ace/evaluate.x b/noao/nproto/ace/evaluate.x new file mode 100644 index 00000000..c3b5b608 --- /dev/null +++ b/noao/nproto/ace/evaluate.x @@ -0,0 +1,641 @@ +include +include +include +include "ace.h" +include "cat.h" +include "objs.h" +include "evaluate.h" + + +# EVALUATE -- Evaluate object parameters. + +procedure evaluate (evl, cat, im, om, skymap, sigmap, gainmap, expmap, logfd) + +pointer evl #I Parameters +pointer cat #I Catalog structure +pointer im #I Image pointer +pointer om #I Object mask pointer +pointer skymap #I Sky map +pointer sigmap #I Sigma map +pointer gainmap #I Gain map +pointer expmap #I Exposure map +int logfd #I Logfile + +int i, n, c, l, nc, nl, c1, c2, nummax, num, nobjsap +real x, x2, y, y2, s, s2, f, f2, val, sky, ssig, s2x, s2y +pointer objs, obj, rlptr +pointer data, skydata, ssigdata, gaindata, expdata, sigdata +pointer sp, v, rl, sum_s2x, sum_s2y + +int andi(), ori(), ctor() +real imgetr() +bool pm_linenotempty() +errchk salloc, calloc, malloc, evgdata + +begin + call smark (sp) + call salloc (v, PM_MAXDIM, TY_LONG) + call salloc (rl, 3+3*IM_LEN(im,1), TY_INT) + + if (logfd != NULL) + call fprintf (logfd, " Evaluate objects:\n") + + objs = CAT_OBJS(cat) + nummax = CAT_NUMMAX(cat) + + nc = IM_LEN(im,1) + nl = IM_LEN(im,2) + + # Allocate work arrays. + call salloc (sigdata, nc, TY_REAL) + call salloc (sum_s2x, nummax, TY_REAL) + call salloc (sum_s2y, nummax, TY_REAL) + call aclrr (Memr[sum_s2x], nummax) + call aclrr (Memr[sum_s2y], nummax) + + # Initialize isophotal quantities. + do i = NUMSTART-1, nummax-1 { + obj = Memi[objs+i] + if (obj == NULL) + next + OBJ_NPIX(obj) = 0 + OBJ_SKY(obj) = 0. + OBJ_PEAK(obj) = 0. + OBJ_FLUX(obj) = 0. + OBJ_X1(obj) = 0. + OBJ_Y1(obj) = 0. + OBJ_X2(obj) = 0. + OBJ_Y2(obj) = 0. + OBJ_XY(obj) = 0. + OBJ_SIG(obj) = 0. + OBJ_ISIGAVG(obj) = 0. + OBJ_ISIGAVG2(obj) = INDEFR + OBJ_FLUXVAR(obj) = 0. + OBJ_XVAR(obj) = 0. + OBJ_YVAR(obj) = 0. + OBJ_XYCOV(obj) = 0. + } + + # Initialize aperture photometry. + call evapinit (cat, nobjsap) + + # Get magnitude zero. + if (EVL_MAGZERO(evl,1) == '!') { + iferr (CAT_MAGZERO(cat) = imgetr (im, EVL_MAGZERO(evl,2))) { + call erract (EA_WARN) + CAT_MAGZERO(cat) = INDEFR + } + } else { + i = 1 + if (ctor (EVL_MAGZERO(evl,1), i, CAT_MAGZERO(cat)) == 0) + CAT_MAGZERO(cat) = INDEFR + } + call catputr (cat, "magzero", CAT_MAGZERO(cat)) + + # Go through the lines of the image accumulating the image data + # into the parameters. The data is read the first time it is + # required. + Memi[v] = 1 + do l = 1, nl { + Memi[v+1] = l + data = NULL + + # Do circular aperture photometry. Check nobjsap to avoid + # subroutine call. + if (nobjsap > 0) + call evapeval (l, im, skymap, sigmap, gainmap, expmap, + data, skydata, ssigdata, gaindata, expdata, sigdata) + + # Accumulate object region quantities if there are object + # regions in the current line. + if (!pm_linenotempty (om, Memi[v])) + next + call pmglri (om, Memi[v], Memi[rl], 0, nc, 0) + + # Go through each object region. + rlptr = rl + do i = 2, Memi[rl] { + rlptr = rlptr + 3 + c1 = Memi[rlptr] + c2 = c1 + Memi[rlptr+1] - 1 + num = MNUM(Memi[rlptr+2]) + + # Do all unevaluated objects and their parents. + while (num >= NUMSTART) { + if (data == NULL) + call evgdata (l, im, skymap, sigmap, gainmap, expmap, + data, skydata, ssigdata, gaindata, expdata, sigdata) + + obj = Memi[objs+num-1] + if (obj == NULL) + break + + if (OBJ_NPIX(obj) == 0) { + val = Memr[data+c1-1] + sky = Memr[skydata+c1-1] + ssig = Memr[ssigdata+c1-1] + + OBJ_XMIN(obj) = c1 + OBJ_XMAX(obj) = c1 + OBJ_YMIN(obj) = l + OBJ_YMAX(obj) = l + OBJ_ISIGMAX(obj) = (val - sky) / ssig + } + + s2x = Memr[sum_s2x+num-1] + s2y = Memr[sum_s2y+num-1] + do c = c1, c2 { + val = Memr[data+c-1] + sky = Memr[skydata+c-1] + ssig = Memr[ssigdata+c-1] + s = Memr[sigdata+c-1] + + x = c - OBJ_XMIN(obj) + y = l - OBJ_YMIN(obj) + x2 = x * x + y2 = y * y + s2 = s * s + + OBJ_NPIX(obj) = OBJ_NPIX(obj) + 1 + OBJ_SKY(obj) = OBJ_SKY(obj) + sky + OBJ_SIG(obj) = OBJ_SIG(obj) + ssig + val = val - sky + if (val > OBJ_PEAK(obj)) + OBJ_PEAK(obj) = val + OBJ_FLUX(obj) = OBJ_FLUX(obj) + val + OBJ_FLUXVAR(obj) = OBJ_FLUXVAR(obj) + s2 + + OBJ_XMIN(obj) = min (OBJ_XMIN(obj), c) + OBJ_XMAX(obj) = max (OBJ_XMAX(obj), c) + OBJ_X1(obj) = OBJ_X1(obj) + x * val + OBJ_X2(obj) = OBJ_X2(obj) + x2 * val + OBJ_XVAR(obj) = OBJ_XVAR(obj) + x2 * s2 + s2x = s2x + x * s2 + + OBJ_YMIN(obj) = min (OBJ_YMIN(obj), l) + OBJ_YMAX(obj) = max (OBJ_YMAX(obj), l) + OBJ_Y1(obj) = OBJ_Y1(obj) + y * val + OBJ_Y2(obj) = OBJ_Y2(obj) + y2 * val + OBJ_YVAR(obj) = OBJ_YVAR(obj) + y2 * s2 + s2y = s2y + y * s2 + + OBJ_XY(obj) = OBJ_XY(obj) + x * y * val + OBJ_XYCOV(obj) = OBJ_XYCOV(obj) + x * y * s2 + + val = val / ssig + OBJ_ISIGAVG(obj) = OBJ_ISIGAVG(obj) + val + OBJ_ISIGMAX(obj) = max (OBJ_ISIGMAX(obj), val) + + } + Memr[sum_s2x+num-1] = s2x + Memr[sum_s2y+num-1] = s2y + + num = OBJ_PNUM(obj) + } + } + } + + # Finish up the evaluations. + do i = NUMSTART-1, nummax-1 { + obj = Memi[objs+i] + if (obj == NULL) + next + n = OBJ_NPIX(obj) + if (n > 0) { + OBJ_SKY(obj) = OBJ_SKY(obj) / n + f = OBJ_FLUX(obj) + if (f > 0.) { + f2 = f * f + x = OBJ_X1(obj) / f + s2x = Memr[sum_s2x+i] + s2y = Memr[sum_s2y+i] + + OBJ_X1(obj) = x + OBJ_XMIN(obj) + OBJ_X2(obj) = OBJ_X2(obj) / f - x * x + OBJ_XVAR(obj) = (OBJ_XVAR(obj) - 2 * x * s2x + + x * x * OBJ_FLUXVAR(obj)) / f2 + + y = OBJ_Y1(obj) / f + OBJ_Y1(obj) = y + OBJ_YMIN(obj) + OBJ_Y2(obj) = OBJ_Y2(obj) / f - y * y + OBJ_YVAR(obj) = (OBJ_YVAR(obj) - 2 * y * s2y + + y * y * OBJ_FLUXVAR(obj)) / f2 + + OBJ_XY(obj) = OBJ_XY(obj) / f - x * y + OBJ_XYCOV(obj) = (OBJ_XYCOV(obj) - x * s2x - + y * s2y + x * y * OBJ_FLUXVAR(obj)) / f2 + + if (IS_INDEFR(OBJ_XAP(obj))) + OBJ_XAP(obj) = OBJ_X1(obj) + if (IS_INDEFR(OBJ_YAP(obj))) + OBJ_YAP(obj) = OBJ_Y1(obj) + } else { + OBJ_X1(obj) = INDEFR + OBJ_Y1(obj) = INDEFR + OBJ_X2(obj) = INDEFR + OBJ_Y2(obj) = INDEFR + OBJ_XY(obj) = INDEFR + OBJ_XVAR(obj) = INDEFR + OBJ_YVAR(obj) = INDEFR + OBJ_XYCOV(obj) = INDEFR + OBJ_FLUXVAR(obj) = INDEFR + } + if (OBJ_PEAK(obj) == 0.) + OBJ_PEAK(obj) = INDEFR + OBJ_SIG(obj) = OBJ_SIG(obj) / n + OBJ_ISIGAVG(obj) = OBJ_ISIGAVG(obj) / sqrt(real(n)) + } + SETFLAG (obj, OBJ_EVAL) + } + + # Do aperture photometry if we had to wait for the aperture centers + # to be defined. + if (nobjsap == 0) { + call evapinit (cat, nobjsap) + if (nobjsap > 0) { + Memi[v] = 1 + do l = 1, nl { + Memi[v+1] = l + data = NULL + call evapeval (l, im, skymap, sigmap, gainmap, expmap, + data, skydata, ssigdata, gaindata, expdata, sigdata) + } + } + } + call evapfree () + + # Set apportioned fluxes. + call evapportion (cat, Memr[sum_s2x]) + + # Set WCS coordinates. + call evalwcs (cat, im) + + call sfree (sp) +end + + +# EVAPINIT -- Initialize aperture photometry. nobjsap will signal whether +# there are any objects to evaluate. + +procedure evapinit (cat, nobjsap) + +pointer cat #I Catalog +int nobjsap #O Number of objects for aperture evaluation + +int i, nummax +pointer tbl, stp, sym, apflux, obj, sthead(), stnext() + +int ycompare() +extern ycompare +errchk calloc, malloc + +int nobjs # Number of objects to evaluate +int naps # Number of apertures per object +real rmax # Maximum aperture radius +pointer r2aps # Array of aperture radii squared (ptr) +pointer ysort # Array of Y sorted object number indices (ptr) +int ystart # Index of first object to consider +pointer objs # Array of object structure (ptr) +common /evapcom/ nobjs, naps, rmax, r2aps, ysort, ystart, objs + +begin + nobjsap = 0 + nobjs = 0 + naps = 0 + r2aps = NULL + ysort = NULL + + tbl = CAT_OUTTBL(cat) + if (tbl == NULL) + return + stp = TBL_STP(tbl) + + # Determine number of apertures. + naps = 0 + for (sym=sthead(stp); sym!=NULL; sym=stnext(stp,sym)) { + if (ENTRY_ID(sym) != ID_APFLUX) + next + } + if (naps == 0) + return + + objs = CAT_OBJS(cat) + nummax = CAT_NUMMAX(cat) + + # Allocate memory. + call calloc (CAT_APFLUX(cat), nummax*naps, TY_REAL) + call malloc (r2aps, naps, TY_REAL) + call malloc (ysort, nummax, TY_INT) + + # Get the maximum radius since that will define the line + # limits needed for each object. Compute array of radius squared + # for the apertures. Pixels are checked for being in the aperture + # in r^2 to avoid square roots. + rmax = 0. + naps = 0 + for (sym=sthead(stp); sym!=NULL; sym=stnext(stp,sym)) { + if (ENTRY_ID(sym) != ID_APFLUX) + next + rmax = max (ENTRY_RAP(sym), rmax) + Memr[r2aps+naps] = ENTRY_RAP(sym) ** 2 + naps = naps + 1 + } + + # Allocate regions of the apflux array to objects with + # defined aperture centers. For the objects create a sorted + # index array by YAP so that we can quickly find objects + # which include a particular line in their apertures. + + apflux = CAT_APFLUX(cat) + do i = NUMSTART-1, nummax-1 { + obj = Memi[objs+i] + if (obj == NULL) + next + if (IS_INDEFR(OBJ_XAP(obj)) || IS_INDEFR(OBJ_YAP(obj))) + next + OBJ_APFLUX(obj) = apflux + apflux = apflux + naps + Memi[ysort+nobjsap] = i + nobjsap = nobjsap + 1 + } + + if (nobjsap > 1) + call gqsort (Memi[ysort], nobjsap, ycompare, objs) + + if (nobjsap == 0) { + call mfree (CAT_APFLUX(cat), TY_REAL) + call evapfree () + } +end + + +# EVAPFREE -- Free aperture photometry memory. + +procedure evapfree () + +int nobjs # Number of objects to evaluate +int naps # Number of apertures per object +real rmax # Maximum aperture radius +pointer r2aps # Array of aperture radii squared (ptr) +pointer ysort # Array of Y sorted object number indices (ptr) +int ystart # Index of first object to consider +pointer objs # Array of object structure (ptr) +common /evapcom/ nobjs, naps, rmax, r2aps, ysort, ystart, objs + +begin + call mfree (r2aps, TY_REAL) + call mfree (ysort, TY_INT) +end + + +# EVAPEVAL -- Do circular aperture photometry. Maintain i1 as the +# first entry in the sorted index array to be considered. All +# earlier entries will have all aperture lines less than the +# current line. Break on the first object whose minimum aperture +# line is greater than the current line. + +procedure evapeval (l, im, skymap, sigmap, gainmap, expmap, data, skydata, + ssigdata, gaindata, expdata, sigdata) + +int l #I Line +pointer im #I Image +pointer skymap #I Sky map +pointer sigmap #I Sigma map +pointer gainmap #I Gain map +pointer expmap #I Exposure map +pointer data #O Image data +pointer skydata #O Sky data +pointer ssigdata #O Sky sigma data +pointer gaindata #O Gain data +pointer expdata #O Exposure data +pointer sigdata #O Total sigma data + +int i, j, nc, c +real x, y, l2, r2, val, sky +pointer obj, apflux + +int nobjs # Number of objects to evaluate +int naps # Number of apertures per object +real rmax # Maximum aperture radius +pointer r2aps # Array of aperture radii squared (ptr) +pointer ysort # Array of Y sorted object number indices (ptr) +int ystart # Index of first object to consider +pointer objs # Array of object structure (ptr) +common /evapcom/ nobjs, naps, rmax, r2aps, ysort, ystart, objs + +begin + nc = IM_LEN(im,1) + do i = ystart, nobjs { + obj = Memi[objs+Memi[ysort+i-1]] + y = OBJ_YAP(obj) + if (y - rmax > l) + break + if (y + rmax < l) { + ystart = ystart + 1 + next + } + x = OBJ_XAP(obj) + apflux = OBJ_APFLUX(obj) + if (data == NULL) + call evgdata (l, im, skymap, sigmap, gainmap, expmap, + data, skydata, ssigdata, gaindata, expdata, sigdata) + + # Accumulate data within in the apertures using the r^2 + # values. Currently partial pixels are not considered and + # errors are not evaluated. + # Note that bad pixels or object overlaps are not excluded + # in the apertures. + l2 = (l - y) ** 2 + do c = max (0, int(x-rmax)), min (nc, int(x+rmax+1)) { + r2 = (c - x) ** 2 + l2 + do j = 0, naps-1 { + if (r2 < Memr[r2aps+j]) { + val = Memr[data+c-1] + sky = Memr[skydata+c-1] + Memr[apflux+j] = Memr[apflux+j] + (val - sky) + } + } + } + } +end + + +# EVAPPORTION -- Compute apportioned fluxes after the object isophotoal +# fluxes have been computed. + +procedure evapportion (cat, sum_flux) + +pointer cat #I Catalog +real sum_flux[ARB] #I Work array of size NUMMAX + +int nummax, num, pnum, nindef +pointer objs, obj, pobj + +begin + objs = CAT_OBJS(cat) + nummax = CAT_NUMMAX(cat) + + call aclrr (sum_flux, nummax) + do num = NUMSTART, nummax { + obj = Memi[objs+num-1] + if (obj == NULL) + next + pnum = OBJ_PNUM(obj) + if (pnum == 0) { + OBJ_FRAC(obj) = 1. + OBJ_FRACFLUX(obj) = OBJ_FLUX(obj) + next + } + + sum_flux[pnum] = sum_flux[pnum] + max (0., OBJ_FLUX(obj)) + OBJ_FRACFLUX(obj) = INDEFR + } + + nindef = 0 + do num = NUMSTART, nummax { + obj = Memi[objs+num-1] + if (obj == NULL) + next + pnum = OBJ_PNUM(obj) + if (pnum == 0) + next + pobj = Memi[objs+pnum-1] + + if (sum_flux[pnum] > 0.) { + OBJ_FRAC(obj) = max (0., OBJ_FLUX(obj)) / sum_flux[pnum] + if (IS_INDEFR(OBJ_FRACFLUX(pobj))) + nindef = nindef + 1 + else + OBJ_FRACFLUX(obj) = OBJ_FRACFLUX(pobj) * OBJ_FRAC(obj) + } else { + OBJ_FRAC(obj) = INDEFR + OBJ_FRACFLUX(obj) = OBJ_FLUX(obj) + } + } + + while (nindef > 0) { + nindef = 0 + do num = NUMSTART, nummax { + obj = Memi[objs+num-1] + if (obj == NULL) + next + pnum = OBJ_PNUM(obj) + if (pnum == 0) + next + + pobj = Memi[objs+pnum-1] + if (IS_INDEFR(OBJ_FRACFLUX(pobj))) + nindef = nindef + 1 + else { + if (IS_INDEFR(OBJ_FRAC(obj))) + OBJ_FRACFLUX(obj) = OBJ_FLUX(obj) + else + OBJ_FRACFLUX(obj) = OBJ_FRACFLUX(pobj) * OBJ_FRAC(obj) + } + } + } +end + + +# EVALWCS -- Set WCS coordinates. + +procedure evalwcs (cat, im) + +pointer cat #I Catalog structure +pointer im #I IMIO pointer + +int i +pointer mw, ct, objs, obj, mw_openim(), mw_sctran() +errchk mw_openim + +begin + mw = mw_openim (im) + ct = mw_sctran (mw, "logical", "world", 03B) + + objs = CAT_OBJS(cat) + do i = NUMSTART-1, CAT_NUMMAX(cat)-1 { + obj = Memi[objs+i] + if (obj == NULL) + next + if (IS_INDEFR(OBJ_XAP(obj)) || IS_INDEFR(OBJ_YAP(obj))) { + OBJ_WX(obj) = INDEFD + OBJ_WY(obj) = INDEFD + } else + call mw_c2trand (ct, double(OBJ_XAP(obj)), + double(OBJ_YAP(obj)), OBJ_WX(obj), OBJ_WY(obj)) + } + + call mw_ctfree (ct) + call mw_close (mw) +end + + +# YCOMPARE -- Compare Y values of two objects for sorting. + +int procedure ycompare (objs, i1, i2) + +pointer objs #I Pointer to array of objects +int i1 #I Index of first object to compare +int i2 #I Index of second object to compare + +real y1, y2 + +begin + y1 = OBJ_YAP(Memi[objs+i1]) + y2 = OBJ_YAP(Memi[objs+i2]) + if (y1 < y2) + return (-1) + else if (y1 > y2) + return (1) + else + return (0) +end + + +# EVGDATA -- Get evaluation data for an image line. + +procedure evgdata (l, im, skymap, sigmap, gainmap, expmap, data, skydata, + ssigdata, gaindata, expdata, sigdata) + +int l #I Line +pointer im #I Image +pointer skymap #I Sky map +pointer sigmap #I Sigma map +pointer gainmap #I Gain map +pointer expmap #I Exposure map +pointer data #O Image data +pointer skydata #O Sky data +pointer ssigdata #O Sky sigma data +pointer gaindata #O Gain data +pointer expdata #O Exposure data +pointer sigdata #O Total sigma data + +int nc +pointer imgl2r(), map_glr() +errchk imgl2r, map_glr, noisemodel + +begin + nc = IM_LEN(im,1) + data = imgl2r (im, l) + skydata = map_glr (skymap, l, READ_ONLY) + ssigdata = map_glr (sigmap, l, READ_ONLY) + if (gainmap == NULL && expmap == NULL) + sigdata = ssigdata + else if (expmap == NULL) { + gaindata = map_glr (gainmap, l, READ_ONLY) + call noisemodel (Memr[data], Memr[skydata], + Memr[ssigdata], Memr[gaindata], INDEFR, + Memr[sigdata], nc) + } else if (gainmap == NULL) { + expdata = map_glr (expmap, l, READ_WRITE) + call noisemodel (Memr[data], Memr[skydata], + Memr[ssigdata], INDEFR, Memr[expdata], + Memr[sigdata], nc) + } else { + gaindata = map_glr (gainmap, l, READ_ONLY) + expdata = map_glr (expmap, l, READ_WRITE) + call noisemodel (Memr[data], Memr[skydata], + Memr[ssigdata], Memr[gaindata], + Memr[expdata], Memr[sigdata], nc) + } +end -- cgit