diff options
Diffstat (limited to 'noao/digiphot/apphot/center')
28 files changed, 3383 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/center/apbcenter.x b/noao/digiphot/apphot/center/apbcenter.x new file mode 100644 index 00000000..9a607ee3 --- /dev/null +++ b/noao/digiphot/apphot/center/apbcenter.x @@ -0,0 +1,96 @@ +include <fset.h> +include "../lib/apphot.h" +include "../lib/display.h" + +# APBCENTER -- Procedure to determine accurate centers for a list of objects +# in batch mode using a simple text file list. + +procedure apbcenter (ap, im, cl, out, id, ld, mgd, gid, interactive) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the IRAF image +int cl # the coordinate file descriptor +int out # output file descriptor +int id # output file sequence number +int ld # the coordinate file list number +pointer mgd # pointer to the metacode file stream +pointer gid # pointer to the image display stream +int interactive # mode of use + +int stdin, ier, ild +pointer sp, str +real wx, wy +int apfitcenter(), fscan(), nscan(), strncmp(), apstati() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call fstats (cl, F_FILENAME, Memc[str], SZ_FNAME) + + # Initialize. + ild = ld + + # Print query. + if (strncmp ("STDIN", Memc[str], 5) == 0) { + stdin = YES + call printf ("Type object x and y coordinates (^D or ^Z to end): ") + call flush (STDOUT) + } else + stdin = NO + + # Loop over the coordinate file. + while (fscan (cl) != EOF) { + + # Get and store the coordinates. + call gargr (wx) + call gargr (wy) + if (nscan () != 2) { + if (stdin == YES) { + call printf ( + "Type object x and y coordinates (^D or ^Z to end): ") + call flush (STDOUT) + } + next + } + + # Transform the input coordinates. + switch (apstati(ap,WCSIN)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_itol (ap, wx, wy, wx, wy, 1) + case WCS_TV: + call ap_vtol (im, wx, wy, wx, wy, 1) + default: + ; + } + + call apsetr (ap, CWX, wx) + call apsetr (ap, CWY, wy) + + # Fit the center. + ier = apfitcenter (ap, im, wx, wy) + + # Print and/or plot the results. + if (interactive == YES) { + call ap_qcenter (ap, ier) + if (gid != NULL) + call apmark (ap, gid, apstati (ap, MKCENTER), NO, NO) + } + call ap_cplot (ap, id, mgd, YES) + if (id == 1) + call ap_param (ap, out, "center") + call ap_pcenter (ap, out, id, ild, ier) + + # Setup for the next object. + id = id + 1 + ild = ild + 1 + call apsetr (ap, WX, wx) + call apsetr (ap, WY, wy) + if (stdin == YES) { + call printf ( + "Type object x and y coordinates (^D or ^Z to end): ") + call flush (STDOUT) + } + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/center/apcconfirm.x b/noao/digiphot/apphot/center/apcconfirm.x new file mode 100644 index 00000000..6c570bd4 --- /dev/null +++ b/noao/digiphot/apphot/center/apcconfirm.x @@ -0,0 +1,75 @@ +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/center.h" + +# AP_CCONFIRM -- Procedure to confirm the critical center parameters. + +procedure ap_cconfirm (ap, out, stid) + +pointer ap # pointer to the apphot structure +int out # output file descriptor +int stid # output file sequence number + +pointer sp, str +real fwhmpsf, capert, skysigma, datamin, datamax +int apstati() +real apstatr(), ap_vfwhmpsf(), ap_vcapert(), ap_vsigma() +real ap_vdatamin(), ap_vdatamax() + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + + call printf ("\n") + + # Confirm the centering algorithm. + call ap_vcstring (ap, Memc[str], SZ_FNAME) + + if (apstati (ap, CENTERFUNCTION) != AP_NONE) { + + # Confirm the fwhmpsf. + if (apstati (ap, CENTERFUNCTION) != AP_CENTROID1D) { + fwhmpsf = ap_vfwhmpsf (ap) + } else + fwhmpsf = apstatr (ap, FWHMPSF) + + # Confirm the centering box. + capert = 2.0 * ap_vcapert (ap) + + # Confirm the sky sigma parameter. + skysigma = ap_vsigma (ap) + + } else { + + fwhmpsf = apstatr (ap, FWHMPSF) + capert = 2.0 * apstatr (ap, CAPERT) + skysigma = apstatr (ap, SKYSIGMA) + + } + + # Confirm the sky sigma parameter. + + # Confirm the minimum and maximum good data values. + datamin = ap_vdatamin (ap) + datamax = ap_vdatamax (ap) + + call printf ("\n") + + # Update the database file. + if (out != NULL && stid > 1) { + call ap_sparam (out, KY_CSTRING, Memc[str], UN_CALGORITHM, + "centering algorithm") + call ap_rparam (out, KY_FWHMPSF, fwhmpsf, UN_ASCALEUNIT, + "full width half maximum of the psf") + call ap_rparam (out, KY_CAPERT, capert, UN_CSCALEUNIT, + "centering box width") + call ap_rparam (out, KY_SKYSIGMA, skysigma, UN_NCOUNTS, + "standard deviation of 1 sky pixel") + call ap_rparam (out, KY_DATAMIN, datamin, UN_ACOUNTS, + "minimum good data value") + call ap_rparam (out, KY_DATAMAX, datamax, UN_ACOUNTS, + "maximum good data value") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/center/apcencolon.x b/noao/digiphot/apphot/center/apcencolon.x new file mode 100644 index 00000000..9f2f40f8 --- /dev/null +++ b/noao/digiphot/apphot/center/apcencolon.x @@ -0,0 +1,329 @@ +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/center.h" + +# APCENTERCOLON -- Process colon commands from the centering task. + +procedure apcentercolon (ap, im, cl, out, stid, ltid, cmdstr, newimage, + newbuf, newfit) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the iraf image +int cl # coordinate file descriptor +int out # output file descriptor +int stid # output file sequence number +int ltid # input list sequence number +char cmdstr[ARB] # command string +int newimage # new image ? +int newbuf # new center buffer ? +int newfit # new center fit ? + +int junk +pointer sp, incmd, outcmd +int strdic() + +begin + # Allocate working space. + call smark (sp) + call salloc (incmd, SZ_LINE, TY_CHAR) + call salloc (outcmd, SZ_LINE, TY_CHAR) + + # Get the command. + call sscan (cmdstr) + call gargwrd (Memc[incmd], SZ_LINE) + if (Memc[incmd] == EOS) { + call sfree (sp) + return + } + + # Process the command. + if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, CCMDS) != 0) + call ap_ccolon (ap, out, stid, cmdstr, newbuf, newfit) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, APCMDS) != 0) + call ap_apcolon (ap, im, cl, out, stid, ltid, cmdstr, newimage, + newbuf, newfit, junk, junk, junk, junk) + else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, NCMDS) != 0) + call ap_nscolon (ap, im, out, stid, cmdstr, newbuf, newfit, + junk, junk, junk, junk) + else + call ap_cimcolon (ap, cmdstr) + + call sfree (sp) +end + + +# AP_CCOLON -- Process colon commands affecting the centering parameters. + +procedure ap_ccolon (ap, out, stid, cmdstr, newbuf, newfit) + +pointer ap # pointer to the apphot structure +int out # output file descriptor +int stid # output list id number +char cmdstr # command string +int newbuf # new sky buffer +int newfit # new sky fit + +bool bval +int ival, ncmd, stat +pointer sp, cmd, str +real rval + +bool itob() +int nscan(), strdic(), btoi(), apstati() +real apstatr() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (str, SZ_FNAME, TY_CHAR) + + # Get the command. + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call sfree (sp) + return + } + + # Process the command. + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, CCMDS) + switch (ncmd) { + + case CCMD_CBOXWIDTH: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_CAPERT) + call pargr (2.0 * apstatr (ap, CAPERT)) + call pargstr (UN_CSCALEUNIT) + } else { + call apsetr (ap, CAPERT, (rval / 2.0)) + if (stid > 1) + call ap_rparam (out, KY_CAPERT, rval, UN_CSCALEUNIT, + "centering subraster width") + newbuf = YES + newfit = YES + } + + case CCMD_CALGORITHM: + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call apstats (ap, CSTRING, Memc[str], SZ_FNAME) + call printf ("%s = %s %s\n") + call pargstr (KY_CSTRING) + call pargstr (Memc[str]) + call pargstr (UN_CALGORITHM) + } else { + stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, CFUNCS) + if (stat > 0) { + call apseti (ap, CENTERFUNCTION, stat) + call apsets (ap, CSTRING, Memc[cmd]) + if (stid > 1) + call ap_sparam (out, KY_CSTRING, Memc[cmd], + UN_CALGORITHM, "centering algorihtm") + newbuf = YES + newfit = YES + } + } + + case CCMD_CTHRESHOLD: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_CTHRESHOLD) + call pargr (apstatr (ap, CTHRESHOLD)) + call pargstr (UN_CSIGMA) + } else { + call apsetr (ap, CTHRESHOLD, rval) + if (stid > 1) + call ap_rparam (out, KY_CTHRESHOLD, rval, UN_CSIGMA, + "threshold for centering") + newbuf = YES + newfit = YES + } + + case CCMD_MAXSHIFT: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_MAXSHIFT) + call pargr (apstatr (ap, MAXSHIFT)) + call pargstr (UN_CSCALEUNIT) + } else { + call apsetr (ap, MAXSHIFT, rval) + if (stid > 1) + call ap_rparam (out, KY_MAXSHIFT, rval, UN_CSCALEUNIT, + "maximum shift") + newfit = YES + } + + case CCMD_MINSNRATIO: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g\n") + call pargstr (KY_MINSNRATIO) + call pargr (apstatr (ap, MINSNRATIO)) + } else { + call apsetr (ap, MINSNRATIO, rval) + if (stid > 1) + call ap_rparam (out, KY_MINSNRATIO, rval, + UN_CNUMBER, "minimum signal-to-noise ratio") + newfit = YES + } + + case CCMD_CMAXITER: + call gargi (ival) + if (nscan () == 1) { + call printf ("%s = %d\n") + call pargstr (KY_CMAXITER) + call pargi (apstati (ap, CMAXITER)) + } else { + call apseti (ap, CMAXITER, ival) + if (stid > 1) + call ap_iparam (out, KY_CMAXITER, ival, UN_CNUMBER, + "maximum number of iterations") + newfit = YES + } + + case CCMD_CLEAN: + call gargb (bval) + if (nscan () == 1) { + call printf ("%s = %b\n") + call pargstr (KY_CLEAN) + call pargb (itob (apstati (ap, CLEAN))) + } else { + call apseti (ap, CLEAN, btoi (bval)) + if (stid > 1) + call ap_bparam (out, KY_CLEAN, bval, UN_CSCALEUNIT, + "apply clean algorithm before centering") + newbuf = YES + newfit = YES + } + + case CCMD_RCLEAN: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_RCLEAN) + call pargr (apstatr (ap, RCLEAN)) + call pargstr (UN_CSCALEUNIT) + } else { + call apsetr (ap, RCLEAN, rval) + if (stid > 1) + call ap_rparam (out, KY_RCLEAN, rval, UN_CSCALEUNIT, + "cleaning radius") + newbuf = YES + newfit = YES + } + + case CCMD_RCLIP: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_RCLIP) + call pargr (apstatr (ap, RCLIP)) + call pargstr (UN_CSCALEUNIT) + } else { + call apsetr (ap, RCLIP, rval) + if (stid > 1) + call ap_rparam (out, KY_RCLIP, rval, UN_CSCALEUNIT, + "clipping radius") + newbuf = YES + newfit = YES + } + + case CCMD_KCLEAN: + call gargr (rval) + if (nscan() == 1) { + call printf ("%s = %g %s\n") + call pargstr (KY_SIGMACLEAN) + call pargr (apstatr (ap, SIGMACLEAN)) + call pargstr (UN_CSIGMA) + } else { + call apsetr (ap, SIGMACLEAN, rval) + if (stid > 1) + call ap_rparam (out, KY_SIGMACLEAN, rval, + UN_CSIGMA, "k-sigma clean rejection criterion") + newbuf = YES + newfit = YES + } + + case CCMD_MKCENTER: + call gargb (bval) + if (nscan() == 1) { + call printf ("%s = %b\n") + call pargstr (KY_MKCENTER) + call pargb (itob (apstati (ap, MKCENTER))) + } else + call apseti (ap, MKCENTER, btoi (bval)) + + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + + call sfree (sp) +end + + +# AP_CIMCOLON -- Process colon commands for the center task that do +# not affect the centering parameters + +procedure ap_cimcolon (ap, cmdstr) + +pointer ap # pointer to the apphot structure +char cmdstr[ARB] # command string + +bool bval +int ncmd +pointer sp, cmd +bool itob() +int strdic(), nscan(), apstati(), btoi() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Get the command. + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + if (Memc[cmd] == EOS) { + call sfree (sp) + return + } + + # Process the command. + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, MISC) + switch (ncmd) { + case ACMD_SHOW: + call gargwrd (Memc[cmd], SZ_LINE) + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, CSHOWARGS) + switch (ncmd) { + case CCMD_DATA: + call printf ("\n") + call ap_nshow (ap) + call printf ("\n") + case CCMD_CENTER: + call printf ("\n") + call ap_cpshow (ap) + call printf ("\n") + default: + call printf ("\n") + call ap_cshow (ap) + call printf ("\n") + } + case ACMD_RADPLOTS: + call gargb (bval) + if (nscan () == 1) { + call printf ("%s = %b\n") + call pargstr (KY_RADPLOTS) + call pargb (itob (apstati (ap, RADPLOTS))) + } else { + call apseti (ap, RADPLOTS, btoi (bval)) + } + default: + call printf ("Unknown or ambiguous colon command\7\n") + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/center/apcenter.x b/noao/digiphot/apphot/center/apcenter.x new file mode 100644 index 00000000..8e1b4c68 --- /dev/null +++ b/noao/digiphot/apphot/center/apcenter.x @@ -0,0 +1,357 @@ +include <ctype.h> +include <gset.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/display.h" +include "../lib/center.h" + +define HELPFILE "apphot$center/center.key" + +# APCENTER -- Procedure to determine accurate centers for a list of objects +# interactively. + +int procedure apcenter (ap, im, cl, gd, mgd, id, out, stid, interactive, cache) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the IRAF image +int cl # coordinate list file descriptor +pointer gd # pointer to the graphics stream +pointer mgd # pointer to the plot metacode stream +pointer id # pointer to the image display stream +int out # output file descriptor +int stid # output file sequence number +int interactive # interactive mode +int cache # cache the input pixels + +real wx, wy, xlist, ylist +pointer sp, cmd +int wcs, key, colonkey, newimage, newobject, newcenter, newlist, ier, oid +int ip, prev_num, req_num, ltid, buf_size, old_size, req_size, memstat + +real apstatr() +int ctoi(), clgcur(), apgscur(), apfitcenter(), aprefitcenter(), apstati() +int apgqverify(), apgtverify(), apnew(), ap_memstat(), sizeof() + +define endswitch_ 99 + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Initialize cursor command. + key = ' ' + Memc[cmd] = EOS + + # Initialize fitting parameters. + newimage = NO + newobject = YES + newcenter = YES + ier = AP_OK + + # Initialize sequencing. + newlist = NO + ltid = 0 + + # Loop over the cursor commands. + while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) != + EOF) { + + # Store the current cursor coordinates. + call ap_vtol (im, wx, wy, wx, wy, 1) + call apsetr (ap, CWX, wx) + call apsetr (ap, CWY, wy) + + # Test to see if the cursor has moved. + if (apnew (ap, wx, wy, xlist, ylist, newlist) == YES) { + newobject = YES + newcenter = YES + } + + # Process the colon commands. + switch (key) { + + # Quit. + case 'q': + if (interactive == YES) { + if (apgqverify ("center", ap, key) == YES) { + call sfree (sp) + return (apgtverify (key)) + } + } else { + call sfree (sp) + return (NO) + } + + # Print error messages. + case 'e': + if (interactive == YES) + call ap_cerrors (ap, ier) + + # Get information on keystroke commands. + case '?': + if ((id != NULL) && (id == gd)) + call gpagefile (id, HELPFILE, "") + else if (interactive == YES) + call pagefile (HELPFILE, "[space=cmhelp,q=quit,?=help]") + + # Rewind the coord list. + case 'r': + if (cl != NULL) { + call seek (cl, BOF) + ltid = 0 + } else if (interactive == YES) + call printf ("No coordinate list\n") + + # Draw a centered stellar radial profile plot. + case 'd': + if (interactive == YES) { + call ap_qrad (ap, im, wx, wy, gd) + newobject = YES + newcenter = YES + } + + # Interactively set the centering parameters. + case 'i': + if (interactive == YES) { + call ap_cradsetup (ap, im, wx, wy, gd, out, stid) + newobject = YES + newcenter = YES + } + + # Verify critical parameters. + case 'v': + call ap_cconfirm (ap, out, stid) + newobject = YES + newcenter = YES + + # Measure the next object. + case 'm', 'n': + + # No coordinate file. + if (cl == NULL) { + if (interactive == YES) + call printf ("No coordinate list\n") + goto endswitch_ + } + + # Need to rewind coordinate file. + prev_num = ltid + req_num = ltid + 1 + if (apgscur (cl, id, xlist, ylist, prev_num, req_num, ltid) == + EOF) { + if (interactive == YES) + call printf ( + "End of coordinate list, use r key to rewind\n") + goto endswitch_ + } + + # Convert coordinates if necessary. + switch (apstati (ap, WCSIN)) { + case WCS_PHYSICAL, WCS_WORLD: + call ap_itol (ap, xlist, ylist, xlist, ylist, 1) + case WCS_TV: + call ap_vtol (im, xlist, ylist, xlist, ylist, 1) + default: + ; + } + + # Move to next object. + newlist = YES + if (key == 'm') { + newobject = YES + newcenter = YES + goto endswitch_ + } + + # Measure next object. + ier = apfitcenter (ap, im, xlist, ylist) + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), NO, NO) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_cplot (ap, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qcenter (ap, ier) + if (stid == 1) + call ap_param (ap, out, "center") + call ap_cplot (ap, stid, mgd, YES) + call ap_pcenter (ap, out, stid, ltid, ier) + stid = stid + 1 + newobject = NO + newcenter = NO + + # Save centering parameters in the pset files. + case 'w': + call ap_pcpars (ap) + + # Process apphot : commands. + case ':': + for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1) + ; + colonkey = Memc[cmd+ip-1] + + # Process the apphot :n command. + switch (colonkey) { + case 'm', 'n': + + # Show/set apphot commands + if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') { + call apcentercolon (ap, im, cl, out, stid, ltid, + Memc[cmd], newimage, newobject, newcenter) + goto endswitch_ + } + + # No coordinate list. + if (cl == NULL) { + if (interactive == YES) + call printf ("No coordinate list\n") + goto endswitch_ + } + + # Get next object from list. + ip = ip + 1 + prev_num = ltid + if (ctoi (Memc[cmd], ip, req_num) <= 0) + req_num = ltid + 1 + + # Fetch the next object from the list. + if (apgscur (cl, id, xlist, ylist, prev_num, req_num, + ltid) == EOF) { + if (interactive == YES) + call printf ( + "End of coordinate list, use r key to rewind\n") + goto endswitch_ + } + + # Convert the coordinates. + switch (apstati (ap, WCSIN)) { + case WCS_PHYSICAL, WCS_WORLD: + call ap_itol (ap, xlist, ylist, xlist, ylist, 1) + case WCS_TV: + call ap_vtol (im, xlist, ylist, xlist, ylist, 1) + default: + ; + } + + # Move to the next object. + newlist = YES + if (colonkey == 'm') { + newobject = YES + newcenter = YES + goto endswitch_ + } + + # Measure next object. + ier = apfitcenter (ap, im, xlist, ylist) + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), NO, NO) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_cplot (ap, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qcenter (ap, ier) + if (stid == 1) + call ap_param (ap, out, "center") + call ap_cplot (ap, stid, mgd, YES) + call ap_pcenter (ap, out, stid, ltid, ier) + stid = stid + 1 + newobject = NO + newcenter = NO + + # Show/set apphot parameters. + default: + call apcentercolon (ap, im, cl, out, stid, ltid, Memc[cmd], + newimage, newobject, newcenter) + } + + if (newimage == YES) { + if ((id != NULL) && (id != gd)) + call ap_gswv (id, Memc[cmd], im, 4) + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + } + + newimage = NO + + # Fit the center, and save the results. + case 'f', ' ': + + # Compute the center. + if (newobject == YES) + ier = apfitcenter (ap, im, wx, wy) + else if (newcenter == YES) + ier = aprefitcenter (ap, im, ier) + if (id != NULL) { + call apmark (ap, id, apstati (ap, MKCENTER), NO, NO) + if (id == gd) + call gflush (id) + else + call gframe (id) + } + call ap_cplot (ap, stid, gd, apstati (ap, RADPLOTS)) + if (interactive == YES) + call ap_qcenter (ap, ier) + newobject = NO; newcenter = NO + + # Write the results to the output file(s). + if (key == ' ') { + if (stid == 1) + call ap_param (ap, out, "center") + if (newlist == YES) + call ap_pcenter (ap, out, stid, ltid, ier) + else + call ap_pcenter (ap, out, stid, 0, ier) + call ap_cplot (ap, stid, mgd, YES) + stid = stid + 1 + } + + # Fit centers for the rest of the list. + case 'l': + + # No coordinate list. + if (cl == NULL) { + if (interactive == YES) + call printf ("No coordinate list\n") + goto endswitch_ + } + + # Center rest of coordinate list. + ltid = ltid + 1 + oid = stid + call apbcenter (ap, im, cl, out, stid, ltid, mgd, id, YES) + ltid = ltid + stid - oid + 1 + if (id != NULL) { + if (id == gd) + call gflush (id) + else + call gframe (id) + } + + # Unknown command. + default: + call printf ("Unknown or ambiguous keystroke command\n") + } + +endswitch_ + # Setup for the next object by setting the default keystroke + # command and storing the old cursor coordinates in the + # centering structure. + + key = ' ' + Memc[cmd] = EOS + call apsetr (ap, WX, apstatr (ap, CWX)) + call apsetr (ap, WY, apstatr (ap, CWY)) + + } + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/center/apcerrors.x b/noao/digiphot/apphot/center/apcerrors.x new file mode 100644 index 00000000..4f8c32bb --- /dev/null +++ b/noao/digiphot/apphot/center/apcerrors.x @@ -0,0 +1,32 @@ +include "../lib/center.h" + +# AP_CERRORS -- Procedure to print the centering task error messages. + +procedure ap_cerrrors (ap, ier) + +pointer ap # apphot pointer (unused) +int ier # centering error code + +begin + switch (ier) { + case AP_CTR_NOAREA: + call printf ("There are no pixels in the centering box.\n") + case AP_CTR_OUTOFBOUNDS: + call printf ( + "The centering region is partially outside the image.\n") + case AP_CTR_NTOO_SMALL: + call printf ("The centering box has too few points.\n") + case AP_CTR_SINGULAR: + call printf ("The centering solution is singular.\n") + case AP_CTR_NOCONVERGE: + call printf ("The centering algorithm does not converge.\n") + case AP_CTR_BADSHIFT: + call printf ("The center shift is large.\n") + case AP_CTR_LOWSNRATIO: + call printf ("The signal to noise ratio is low.\n") + case AP_CTR_BADDATA: + call printf ("Bad data in the centering subraster.\n") + default: + call printf ("") + } +end diff --git a/noao/digiphot/apphot/center/apcfree.x b/noao/digiphot/apphot/center/apcfree.x new file mode 100644 index 00000000..7fb219c7 --- /dev/null +++ b/noao/digiphot/apphot/center/apcfree.x @@ -0,0 +1,46 @@ +include "../lib/apphotdef.h" +include "../lib/centerdef.h" + +# APCFREE -- Procedure to free the centering data structure. + +procedure apcfree (ap) + +pointer ap # pointer to the apphot structure + +begin + if (ap == NULL) + return + if (AP_NOISE(ap) != NULL) + call ap_noisecls (ap) + if (AP_PCENTER(ap) != NULL) + call ap_ctrcls (ap) + if (AP_PDISPLAY(ap) != NULL) + call ap_dispcls (ap) + if (AP_IMBUF(ap) != NULL) + call mfree (AP_IMBUF(ap), TY_REAL) + if (AP_MW(ap) != NULL) + call mw_close (AP_MW(ap)) + call mfree (ap, TY_STRUCT) +end + + +# AP_CTRCLS -- Procedure to close up the centering structure arrays. + +procedure ap_ctrcls (ap) + +pointer ap # pointer to apphot structure + +pointer ctr + +begin + ctr = AP_PCENTER(ap) + if (ctr == NULL) + return + if (AP_CTRPIX(ctr) != NULL) + call mfree (AP_CTRPIX(ctr), TY_REAL) + if (AP_XCTRPIX(ctr) != NULL) + call mfree (AP_XCTRPIX(ctr), TY_REAL) + if (AP_YCTRPIX(ctr) != NULL) + call mfree (AP_YCTRPIX(ctr), TY_REAL) + call mfree (AP_PCENTER(ap), TY_STRUCT) +end diff --git a/noao/digiphot/apphot/center/apcinit.x b/noao/digiphot/apphot/center/apcinit.x new file mode 100644 index 00000000..a9f38fd7 --- /dev/null +++ b/noao/digiphot/apphot/center/apcinit.x @@ -0,0 +1,94 @@ +include "../lib/apphotdef.h" +include "../lib/centerdef.h" +include "../lib/center.h" + +# APCINIT - Procedure to initialize the centering structure. + +procedure apcinit (ap, function, cbox, fwhmpsf, noise) + +pointer ap # pointer to the apphot structure +int function # centering algorithm +real cbox # half width of centering box +real fwhmpsf # FWHM of the PSF +int noise # noise function + +begin + # Allocate space. + call malloc (ap, LEN_APSTRUCT, TY_STRUCT) + + # Set up the global apphot package parameters. + call ap_defsetup (ap, fwhmpsf) + + # Setup the noise structure. + call ap_noisesetup (ap, noise) + + # Initialize the centering parameters. + call ap_ctrsetup (ap, function, cbox) + + # Display the options. + call ap_dispsetup (ap) + + # Unused structures are set to null. + AP_PSKY(ap) = NULL + AP_PPSF(ap) = NULL + AP_PPHOT(ap) = NULL + AP_POLY(ap) = NULL + AP_RPROF(ap) = NULL +end + + +# AP_CTRSETUP -- Procedure to setup the centering array parameters. + +procedure ap_ctrsetup (ap, function, cbox) + +pointer ap # pointer to apphot structure +int function # centering function +real cbox # radius of centering aperture + +pointer ctr + +begin + # Allocate space for the centering structure. + call malloc (AP_PCENTER(ap), LEN_CENSTRUCT, TY_STRUCT) + ctr = AP_PCENTER(ap) + AP_CXCUR(ctr) = INDEFR + AP_CYCUR(ctr) = INDEFR + + AP_CENTERFUNCTION(ctr) = function + switch (function) { + case AP_CENTROID1D: + call strcpy ("centroid", AP_CSTRING(ctr), SZ_FNAME) + case AP_GAUSS1D: + call strcpy ("gauss", AP_CSTRING(ctr), SZ_FNAME) + case AP_NONE: + call strcpy ("none", AP_CSTRING(ctr), SZ_FNAME) + case AP_OFILT1D: + call strcpy ("ofilter", AP_CSTRING(ctr), SZ_FNAME) + default: + AP_CENTERFUNCTION(ctr) = DEF_CENTERFUNCTION + call strcpy ("centroid", AP_CSTRING(ctr), SZ_FNAME) + } + AP_CAPERT(ctr) = cbox + AP_CTHRESHOLD(ctr) = DEF_CTHRESHOLD + AP_MINSNRATIO(ctr) = DEF_MINSNRATIO + AP_CMAXITER(ctr) = DEF_CMAXITER + AP_MAXSHIFT(ctr) = DEF_MAXSHIFT + AP_CLEAN(ctr) = DEF_CLEAN + AP_RCLEAN(ctr) = DEF_RCLEAN + AP_RCLIP(ctr) = DEF_RCLIP + AP_SIGMACLEAN(ctr) = DEF_CLEANSIGMA + + AP_NCTRPIX(ctr) = 0 + AP_CTRPIX(ctr) = NULL + AP_XCTRPIX(ctr) = NULL + AP_YCTRPIX(ctr) = NULL + + AP_OXINIT(ctr) = INDEFR + AP_OYINIT(ctr) = INDEFR + AP_XCENTER(ctr) = INDEFR + AP_YCENTER(ctr) = INDEFR + AP_OXCENTER(ctr) = INDEFR + AP_OYCENTER(ctr) = INDEFR + AP_XERR(ctr) = INDEFR + AP_YERR(ctr) = INDEFR +end diff --git a/noao/digiphot/apphot/center/apclean.x b/noao/digiphot/apphot/center/apclean.x new file mode 100644 index 00000000..2e3a827a --- /dev/null +++ b/noao/digiphot/apphot/center/apclean.x @@ -0,0 +1,152 @@ +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/center.h" + +# AP_CLEAN -- Procedure to clean a data array. + +procedure ap_clean (ap, pix, nx, ny, xc, yc) + +pointer ap # pointer to the apphot structure +real pix[nx,ny] # data array +int nx, ny # size of subarray +real xc, yc # center of subarray + +int apstati() +real apstatr() + +begin + switch (apstati (ap, NOISEFUNCTION)) { + case AP_NCONSTANT: + call ap_cclean (pix, nx, ny, xc, yc, apstatr (ap, SCALE) * + apstatr (ap, RCLIP), apstatr (ap, SIGMACLEAN), + apstatr (ap, SKYSIGMA), apstatr (ap, MAXSHIFT)) + case AP_NPOISSON: + call ap_pclean (pix, nx, ny, xc, yc, apstatr (ap, SCALE) * + apstatr (ap, RCLEAN), apstatr (ap, SCALE) * apstatr (ap, + RCLIP), apstatr (ap, SIGMACLEAN), apstatr (ap, SKYSIGMA), + apstatr (ap, EPADU), apstatr (ap, MAXSHIFT)) + default: + return + } +end + + +# AP_CCLEAN -- Procedure to clean the subraster assuming the noise is +# due to a constant sky sigma. + +procedure ap_cclean (pix, nx, ny, cxc, cyc, rclip, kclean, skysigma, + maxshift) + +real pix[nx, ny] # array of pixels +int nx, ny # dimensions of the subarray +real cxc, cyc # initial center +real rclip # cleaning and clipping radius +real kclean # k-sigma clipping factor +real skysigma # maxshift +real maxshift # maximum shift + +int i, ii, ixc, j, jj, jyc, ijunk, jjunk +real mindat, maxdat, rclip2, r2, ksigma + +begin + # Return if indef valued sigma or sigma <= 0. + if (IS_INDEFR(skysigma) || (skysigma <= 0.0)) + return + + # Find the maximum pixel in the subarray and treat this point as + # the center of symmetry if it is less than maxshift from the + # initial center. + + call ap_2dalimr (pix, nx, ny, mindat, maxdat, ijunk, jjunk, ixc, jyc) + if (abs (cxc - ixc) > maxshift || abs (cyc - jyc) > maxshift) { + ixc = int (cxc) + jyc = int (cyc) + } + + # Clip. + rclip2 = rclip ** 2 + ksigma = kclean * skysigma + do j = 1, ny { + jj = 2 * jyc - j + if (jj < 1 || jj > ny) + next + do i = 1, ixc { + ii = 2 * ixc - i + if (ii < 1 || ii > nx) + next + r2 = (i - ixc) ** 2 + (j - jyc) ** 2 + if (r2 > rclip2) { + if (pix[ii,jj] > pix[i,j] + ksigma) + pix[ii,jj] = pix[i,j] + else if (pix[i,j] > pix[ii,jj] + ksigma) + pix[i,j] = pix[ii,jj] + } + } + } +end + + +# AP_PCLEAN -- Procedure to clean the subraster assuming that the noise is +# due to a constant sky value plus poisson noise. + +procedure ap_pclean (pix, nx, ny, cxc, cyc, rclean, rclip, kclean, skysigma, + padu, maxshift) + +real pix[nx, ny] # array of pixels +int nx, ny # dimensions of the subarray +real cxc, cyc # initial center +real rclean, rclip # cleaning and clipping radius +real kclean # k-sigma clipping factor +real skysigma # maxshift +real padu # photons per adu +real maxshift # maximum shift + +int i, ii, ixc, j, jj, jyc, ijunk, jjunk +real mindat, maxdat, rclean2, rclip2, r2, ksigma, ksigma2 + +begin + # Return if indef-valued sigma. + if (IS_INDEFR(skysigma)) + return + + # Find the maximum pixel in the subarray and treat this point as + # the center of symmetry if it is less than maxshift from the + # initial center. + + call ap_2dalimr (pix, nx, ny, mindat, maxdat, ijunk, jjunk, ixc, jyc) + if (abs (cxc - ixc) > maxshift || abs (cyc - jyc) > maxshift) { + ixc = int (cxc) + jyc = int (cyc) + } + + # Clip. + rclean2 = rclean ** 2 + rclip2 = rclip ** 2 + ksigma = kclean * skysigma + ksigma2 = ksigma ** 2 + do j = 1, ny { + jj = 2 * jyc - j + if (jj < 1 || jj > ny) + next + do i = 1, ixc { + ii = 2 * ixc - i + if (ii < 1 || ii > nx) + next + r2 = (i - ixc) ** 2 + (j - jyc) ** 2 + if (r2 > rclean2 && r2 <= rclip2) { + if (pix[ii,jj] > (pix[i,j] + sqrt (pix[i,j] / padu + + ksigma2))) + pix[ii,jj] = pix[i,j] + else if (pix[i,j] > (pix[ii,jj] + sqrt (pix[ii,jj] / padu + + ksigma2))) + pix[i,j] = pix[ii,jj] + } + if (r2 > rclip2) { + if (pix[ii,jj] > pix[i,j] + ksigma) + pix[ii,jj] = pix[i,j] + else if (pix[i,j] > pix[ii,jj] + ksigma) + pix[i,j] = pix[ii,jj] + } + } + } +end diff --git a/noao/digiphot/apphot/center/apcplot.x b/noao/digiphot/apphot/center/apcplot.x new file mode 100644 index 00000000..6e2845c3 --- /dev/null +++ b/noao/digiphot/apphot/center/apcplot.x @@ -0,0 +1,304 @@ +include <mach.h> +include <gset.h> +include <pkg/gtools.h> +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/noise.h" +include "../lib/centerdef.h" +include "../lib/center.h" + +# AP_CPLOT -- Procedure to compute radial profile plots for the centering +# routine. + +procedure ap_cplot (ap, sid, gd, makeplot) + +pointer ap # the pointer to the apphot structure +int sid # the output sequence number of the star +pointer gd # the graphics stream +int makeplot # the make a plot ? + +int nx, ny +pointer ctr, sp, str, r, gt +real xcenter, ycenter, xc, yc, rmin, rmax, imin, imax +real u1, u2, v1, v2, x1, x2, y1, y2 +int apstati() +pointer ap_gtinit() +real apstatr() + +begin + # Check for enabled graphics stream. + if (gd == NULL || makeplot == NO) + return + + # Check for defined center. + xcenter = apstatr (ap, XCENTER) + ycenter = apstatr (ap, YCENTER) + if (IS_INDEFR(xcenter) || IS_INDEFR(ycenter)) + return + + # Check for a centering algorithm. + if (apstati (ap, CENTERFUNCTION) == AP_NONE) + return + + # Get the pixel buffer parameters. + ctr = AP_PCENTER(ap) + nx = AP_CNX(ctr) + ny = AP_CNY(ctr) + if (AP_CTRPIX(ctr) == NULL || nx <= 0 || ny <= 0) + return + xc = AP_CXC(ctr) + (xcenter - apstatr (ap, CXCUR)) + yc = AP_CYC(ctr) + (ycenter - apstatr (ap, CYCUR)) + + # Allocate working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (r, nx * ny, TY_REAL) + + # Compute the radii and the plot limits. + call ap_ijtor2 (Memr[r], nx, ny, xc, yc) + call alimr (Memr[r], nx * ny, rmin, rmax) + call alimr (Memr[AP_CTRPIX(ctr)], nx * ny, imin, imax) + + # Reactivate the work station. + call greactivate (gd, 0) + + # Save old viewport and window coordinates. + call ggview (gd, u1, u2, v1, v2) + call ggwind (gd, x1, x2, y1, y2) + + # Initialize the plot. + call apstats (ap, IMROOT, Memc[str], SZ_LINE) + call sprintf (Memc[str], SZ_LINE, "%s Star %d") + call pargstr (Memc[str]) + call pargi (sid) + gt = ap_gtinit (Memc[str], apstatr (ap,OXINIT), apstatr (ap, OYINIT)) + + # Make the plot. + call gclear (gd) + call ap_cpset (gd, gt, ap, rmin, rmax, imin, imax) + if (apstati (ap, POSITIVE) == YES) + call ap_plotpts (gd, gt, Memr[r], Memr[AP_CTRPIX(ctr)], nx * ny, + rmin, rmax, imin + EPSILONR, imax, "plus") + else + call ap_plotpts (gd, gt, Memr[r], Memr[AP_CTRPIX(ctr)], nx * ny, + rmin, rmax, imin, imax - EPSILONR, "plus") + call ap_cpreset (gd, gt, ap, rmin, rmax, imin, imax) + call ap_cpannotate (gd, ap) + + # Restore the viewport and window coordinates. + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + # Free the space. + call ap_gtfree (gt) + call gdeactivate (gd, 0) + call sfree (sp) +end + + +# AP_CPSET -- Procedure to set up the parameters for the center radial profile +# plot. + +procedure ap_cpset (gd, gt, ap, xmin, xmax, ymin, ymax) + +pointer gd # the graphics stream +pointer gt # the gtools pointer +pointer ap # the apphot pointer +real xmin, xmax # the minimum and maximum radial distance +real ymin, ymax # the minimum and maximum of y axis (ymin not used) + +int fd +pointer sp, str, title +real scale, aspect, datalimit, skysigma, threshold, vx1, vx2, vy1, vy2 +real cthreshold +int stropen(), apstati() +real apstatr(), gstatr() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (title, SZ_LINE, TY_CHAR) + + # Encode the parameter string. + fd = stropen (Memc[str], SZ_LINE, WRITE_ONLY) + call sysid (Memc[title], SZ_LINE) + call fprintf (fd, "%s\n") + call pargstr (Memc[title]) + call fprintf (fd, + "Center: xc=%0.2f yc=%0.2f xerr=%0.2f yerr=%0.2f\n") + call pargr (apstatr (ap, OXCENTER)) + call pargr (apstatr (ap, OYCENTER)) + call pargr (apstatr (ap, XERR)) + call pargr (apstatr (ap, YERR)) + call gt_gets (gt, GTTITLE, Memc[title], SZ_LINE) + call fprintf (fd, "%s\n") + call pargstr (Memc[title]) + call strclose (fd) + + # Store some default plotting parameters. + scale = apstatr (ap, SCALE) + aspect = gstatr (gd, G_ASPECT) + call gsetr (gd, G_ASPECT, 0.75) + + # Set the labels and window. + datalimit = apstatr (ap, CDATALIMIT) + skysigma = apstatr (ap, SKYSIGMA) + cthreshold = apstatr (ap, CTHRESHOLD) + if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold)) + threshold = 0.0 + else + threshold = cthreshold * skysigma + if (apstati (ap, POSITIVE) == YES) { + call gseti (gd, G_XDRAWAXES, 2) + call gswind (gd, xmin / scale, xmax / scale, datalimit + threshold, + datalimit + threshold + ymax) + call glabax (gd, Memc[str], "", "Intensity") + call gseti (gd, G_YDRAWAXES, 0) + call gseti (gd, G_XDRAWAXES, 1) + call ggview (gd, vx1, vx2, vy1, vy2) + call gswind (gd, vx1, vx2, vy1, vy2) + call gswind (gd, xmin, xmax, datalimit + threshold, + datalimit + threshold + ymax) + call glabax (gd, "", + "Radial Distance (lower-pixels, upper-scale units)", "") + } else { + call gseti (gd, G_XDRAWAXES, 2) + call gswind (gd, xmin / scale, xmax / scale, datalimit - + threshold - ymax, datalimit - threshold) + call glabax (gd, Memc[str], "", "Intensity") + call gseti (gd, G_YDRAWAXES, 0) + call gseti (gd, G_XDRAWAXES, 1) + call ggview (gd, vx1, vx2, vy1, vy2) + call gswind (gd, vx1, vx2, vy1, vy2) + call gswind (gd, xmin, xmax, datalimit - threshold - ymax, + datalimit - threshold) + call glabax (gd, "", + "Radial Distance (lower-pixels, upper-scale units)", "") + } + + # Set the window up for plotting the data. + call gseti (gd, G_YDRAWAXES, 3) + call gseti (gd, G_XDRAWAXES, 3) + call gsetr (gd, G_ASPECT, aspect) + call gt_sets (gt, GTTYPE, "mark") + call gt_setr (gt, GTXMIN, xmin) + call gt_setr (gt, GTXMAX, xmax) + if (apstati (ap, POSITIVE) == YES) { + call gt_setr (gt, GTYMIN, ymin) + call gt_setr (gt, GTYMAX, ymax) + } else { + call gt_setr (gt, GTYMIN, ymin) + call gt_setr (gt, GTYMAX, ymax) + } + call gt_swind (gd, gt) + + call sfree (sp) +end + + +# AP_CPANNOTATE -- Procedure to annotate the radial plot in center. + +procedure ap_cpannotate (gd, ap) + +pointer gd # the graphics stream +pointer ap # the apphot structure + +pointer sp, str +real fwhmpsf, capert, datalimit, threshold, skysigma, cthreshold +real xmin, xmax, ymin, ymax +int apstati() +real apstatr() + +begin + # Allocate working space. + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call gseti (gd, G_PLTYPE, GL_DASHED) + call ggwind (gd, xmin, xmax, ymin, ymax) + + fwhmpsf = 0.5 * apstatr (ap, FWHMPSF) * apstatr (ap, SCALE) + capert = 2.0 * fwhmpsf * apstatr (ap, CAPERT) + datalimit = apstatr (ap, CDATALIMIT) + skysigma = apstatr (ap, SKYSIGMA) + cthreshold = apstatr (ap, CTHRESHOLD) + if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold)) + threshold = 0.0 + else + threshold = cthreshold * skysigma + if (apstati (ap, POSITIVE) == YES) + threshold = datalimit + threshold + else + threshold = datalimit - threshold + + # Plot the full width half maximum of the radial profile. + if (fwhmpsf >= xmin && fwhmpsf <= xmax) { + call gamove (gd, fwhmpsf, ymin) + call gadraw (gd, fwhmpsf, ymax) + call sprintf (Memc[str], SZ_LINE, "hwhm = %0.2f") + call pargr (fwhmpsf) + call gtext (gd, fwhmpsf, ymax, Memc[str], "q=h;u=180;v=t;p=r") + } + + # Mark the centering aperture. + if (capert >= xmin && capert <= xmax) { + call gamove (gd, capert, ymin) + call gadraw (gd, capert, ymax) + call sprintf (Memc[str], SZ_LINE, "cbox half-width = %0.2f") + call pargr (capert) + call gtext (gd, capert, ymax, Memc[str], "q=h;u=180;v=t;p=r") + } + + # Mark the threshold level for centering. + call sprintf (Memc[str], SZ_LINE, "threshold = %g") + call pargr (threshold) + call gtext (gd, xmin, ymin, Memc[str], "q=h") + + # Mark the sky sigma if defined. + if (! IS_INDEFR(skysigma) && skysigma >= ymin && skysigma <= ymax) { + call gmark (gd, (xmin + xmax) / 2.0, (ymin + ymax) / 2.0, + GM_VEBAR, -0.25, -skysigma) + call sprintf (Memc[str], SZ_LINE, "sigma = %g") + call pargr (skysigma) + call gtext (gd, (xmin + xmax) / 2.0, (ymin + ymax + skysigma) / 2.0, + Memc[str], "q=h;h=c") + } + + call sfree (sp) +end + + +# AP_CPRSET -- Procedure to reset the plot window after the data points +# have been plotted. + +procedure ap_cpreset (gd, gt, ap, xmin, xmax, ymin, ymax) + +pointer gd # the graphics stream +pointer gt # the gtools pointer +pointer ap # the apphot pointer +real xmin, xmax # the minimum and maximum radial distance +real ymin, ymax # the minimum and maximum of the y axis (ymin not used) + +real datalimit, skysigma, threshold, cthreshold +int apstati() +real apstatr() + +begin + # Set the data window. + datalimit = apstatr (ap, CDATALIMIT) + skysigma = apstatr (ap, SKYSIGMA) + cthreshold = apstatr (ap, CTHRESHOLD) + if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold)) + threshold = 0.0 + else + threshold = cthreshold * skysigma + call gt_setr (gt, GTXMIN, xmin) + call gt_setr (gt, GTXMAX, xmax) + if (apstati (ap, POSITIVE) == YES) { + call gt_setr (gt, GTYMIN, datalimit + threshold) + call gt_setr (gt, GTYMAX, ymax + datalimit + threshold) + } else { + call gt_setr (gt, GTYMIN, datalimit - ymax - threshold) + call gt_setr (gt, GTYMAX, datalimit - threshold) + } + call gt_swind (gd, gt) +end diff --git a/noao/digiphot/apphot/center/apcpshow.x b/noao/digiphot/apphot/center/apcpshow.x new file mode 100644 index 00000000..79fdf43e --- /dev/null +++ b/noao/digiphot/apphot/center/apcpshow.x @@ -0,0 +1,66 @@ +include "../lib/center.h" +include "../lib/display.h" + +# AP_CPSHOW -- Procedure to display the current centering algorithm +# parameters. + +procedure ap_cpshow (ap) + +pointer ap # pointer to the apphot strucuture + +pointer sp, str +bool itob() +int apstati() +real apstatr() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Print the centering algorithm. + call printf ("Centering Parameters\n") + call apstats (ap, CSTRING, Memc[str], SZ_FNAME) + call printf (" %s = %s %s\n") + call pargstr (KY_CSTRING) + call pargstr (Memc[str]) + call pargstr (UN_CALGORITHM) + + # Print the rest of the centering parameters. + call printf (" %s = %g %s %s = %g %s %s = %g\n") + call pargstr (KY_CAPERT) + call pargr (2.0 * apstatr (ap, CAPERT)) + call pargstr (UN_CSCALEUNIT) + call pargstr (KY_CTHRESHOLD) + call pargr (apstatr (ap, CTHRESHOLD)) + call pargstr (UN_CSIGMA) + call pargstr (KY_MINSNRATIO) + call pargr (apstatr (ap, MINSNRATIO)) + + call printf (" %s = %d %s = %g %s\n") + call pargstr (KY_CMAXITER) + call pargi (apstati (ap, CMAXITER)) + call pargstr (KY_MAXSHIFT) + call pargr (apstatr (ap, MAXSHIFT)) + call pargstr (UN_CNUMBER) + + call printf (" %s = %b %s = %g %s\n") + call pargstr (KY_CLEAN) + call pargb (itob (apstati (ap, CLEAN))) + call pargstr (KY_SIGMACLEAN) + call pargr (apstatr (ap, SIGMACLEAN)) + call pargstr (UN_CSIGMA) + + call printf (" %s = %g %s %s = %g %s\n") + call pargstr (KY_RCLEAN) + call pargr (apstatr (ap, RCLEAN)) + call pargstr (UN_CSCALEUNIT) + call pargstr (KY_RCLIP) + call pargr (apstatr (ap, RCLIP)) + call pargstr (UN_CSCALEUNIT) + + call printf (" %s = %b\n") + call pargstr (KY_MKCENTER) + call pargb (itob (apstati (ap, MKCENTER))) + + call sfree (sp) +end diff --git a/noao/digiphot/apphot/center/apcradsetup.x b/noao/digiphot/apphot/center/apcradsetup.x new file mode 100644 index 00000000..8902c5ef --- /dev/null +++ b/noao/digiphot/apphot/center/apcradsetup.x @@ -0,0 +1,99 @@ +include "../lib/display.h" + +define HELPFILE "apphot$center/icenter.key" + +# AP_CRADSETUP -- Procedure to set up apphot interactively + +procedure ap_cradsetup (ap, im, wx, wy, gd, out, stid) + +pointer ap # pointer to apphot structure +pointer im # pointero to the IRAF image +real wx, wy # cursor coordinates +pointer gd # pointer to graphics stream +int out # output file descriptor +int stid # output file sequence number + +int cier, wcs, key +pointer sp, cmd +real xcenter, ycenter, xc, yc, rmin, rmax, imin, imax, rval +real u1, u2, v1, v2, x1, x2, y1, y2 + +int ap_showplot(), apfitcenter(), apstati(), clgcur() +real ap_cfwhmpsf(), ap_ccapert(), ap_csigma() +real ap_cdatamin(), ap_cdatamax(), ap_crclean(), ap_crclip() + +begin + # Check for open graphics stream. + if (gd == NULL) + return + call greactivate (gd, 0) + + # Store old viewport and window coordinates. + call ggview (gd, u1, u2, v1, v2) + call ggwind (gd, x1, x2, y1, y2) + + # Make the plot. + if (ap_showplot (ap, im, wx, wy, gd, xcenter, ycenter, rmin, rmax, + imin, imax) == ERR) { + call gdeactivate (gd, 0) + return + } + + # Allocate temporary space. + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + call printf ( + "Waiting for setup menu command (?=help, v=default setup, q=quit):\n") + + while (clgcur ("gcommands", xc, yc, wcs, key, Memc[cmd], + SZ_LINE) != EOF) { + + switch (key) { + case 'q': + break + case '?': + call gpagefile (gd, HELPFILE, "") + case 'f': + rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax) + case 's': + rval = ap_csigma (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'l': + rval = ap_cdatamin (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'u': + rval = ap_cdatamax (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'c': + rval = ap_ccapert (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'n': + rval = ap_crclean (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'p': + rval = ap_crclip (ap, gd, out, stid, rmin, rmax, imin, imax) + case 'v': + rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax) + rval = ap_ccapert (ap, gd, out, stid, rmin, rmax, imin, imax) + rval = ap_csigma (ap, gd, out, stid, rmin, rmax, imin, imax) + #rval = ap_ccthresh (ap, gd, out, stid, rmin, rmax, imin, imax) + default: + call printf ("Unknown or ambiguous keystroke command\007\n") + } + call printf ( + "Waiting for setup menu command (?=help, v=default setup, q=quit):\n") + } + + # Set up interactively. + call printf ( + "Interactive setup is complete [Type w to save parameters].\n") + + # Restore old window and viewport coordinates. + call gsview (gd, u1, u2, v1, v2) + call gswind (gd, x1, x2, y1, y2) + + # Free plots and working space. + call sfree (sp) + call gdeactivate (gd, 0) + + # Fit the center and plot the results. + cier = apfitcenter (ap, im, xcenter, ycenter) + call ap_cplot (ap, 0, gd, apstati (ap, RADPLOTS)) + call ap_qcenter (ap, cier) +end diff --git a/noao/digiphot/apphot/center/apcshow.x b/noao/digiphot/apphot/center/apcshow.x new file mode 100644 index 00000000..8cedbfc2 --- /dev/null +++ b/noao/digiphot/apphot/center/apcshow.x @@ -0,0 +1,19 @@ +include "../lib/display.h" + +# AP_CSHOW -- Procedure to display the current centering parameters. + +procedure ap_cshow (ap) + +pointer ap # pointer to the apphot strucuture + +bool itob() +int apstati() + +begin + call ap_nshow (ap) + call printf ("\n") + call ap_cpshow (ap) + call printf (" %s = %b\n") + call pargstr (KY_RADPLOTS) + call pargb (itob (apstati (ap, RADPLOTS))) +end diff --git a/noao/digiphot/apphot/center/apcsnratio.x b/noao/digiphot/apphot/center/apcsnratio.x new file mode 100644 index 00000000..0fd0e20d --- /dev/null +++ b/noao/digiphot/apphot/center/apcsnratio.x @@ -0,0 +1,92 @@ +include <mach.h> +include "../lib/center.h" +include "../lib/noise.h" + +# AP_CSNRATIO -- Procedure to estimate the signal to noise ratio for centering. + +real procedure ap_csnratio (ap, array, nx, ny, threshold) + +pointer ap # apphot structure +real array[nx,ny] # data array +int nx, ny # subraster dimensions +real threshold # threshold value for snr computation + +int apstati() +real apstatr(), ap_cratio(), ap_pratio() + +begin + switch (apstati (ap, NOISEFUNCTION)) { + case AP_NCONSTANT: + return (ap_cratio (array, nx, ny, threshold, apstatr (ap, + SKYSIGMA))) + case AP_NPOISSON: + return (ap_pratio (array, nx, ny, threshold, apstatr (ap, + SKYSIGMA), apstatr (ap, EPADU))) + default: + return (MAX_REAL) + } +end + + +# AP_CRATIO -- Procedure to estimate the signal to noise ratio in the +# centering aperture assuming that the noise is due to a constant sky +# sigma. These computations are approximate only. + +real procedure ap_cratio (array, nx, ny, sky, noise) + +real array[nx,ny] # centering buffer +int nx, ny # dimensions of the centering buffer +real sky # mean sky value of the data in ADU +real noise # estimate of sky noise in ADU + +int npts +real signal, tnoise +real asumr() + +begin + npts = nx * ny + signal = asumr (array, npts) - npts * sky + if (IS_INDEFR(noise)) + tnoise = 0.0 + else + tnoise = sqrt (npts * noise ** 2) + if (signal <= 0.0) + return (0.0) + else if (tnoise <= 0.0) + return (MAX_REAL) + else + return (signal / tnoise) +end + + +# AP_PRATIO -- Procedure to estimate the signal to noise ratio in the +# centering aperture assuming the noise is due to a constant sky sigma +# and poisson statistics in the image. These computations are +# approximate only. + +real procedure ap_pratio (array, nx, ny, sky, noise, padu) + +real array[nx,ny] # centering buffer +int nx, ny # dimensions of the centering buffer +real sky # mean sky value of the data in ADU +real noise # estimate of sky noise in ADU +real padu # photons per adu + +int npts +real signal, tnoise +real asumr() + +begin + npts = nx * ny + signal = asumr (array, npts) - npts * sky + if (IS_INDEFR(noise)) + tnoise = sqrt (abs (signal / padu)) + else + tnoise = sqrt (abs (signal / padu) + npts * noise ** 2) + if (signal <= 0.0) + return (0.0) + else if (tnoise <= 0.0) + return (MAX_REAL) + else + return (signal / tnoise) +end diff --git a/noao/digiphot/apphot/center/apctr1d.x b/noao/digiphot/apphot/center/apctr1d.x new file mode 100644 index 00000000..10ee109e --- /dev/null +++ b/noao/digiphot/apphot/center/apctr1d.x @@ -0,0 +1,113 @@ +include <mach.h> +include "../lib/center.h" + +# AP_CTR1D -- Procedure to compute the center from the 1D marginal +# distributions. + +int procedure ap_ctr1d (ctrpix, nx, ny, norm, xc, yc, xerr, yerr) + +real ctrpix[nx, ny] # object to be centered +int nx, ny # dimensions of subarray +real norm # the normalization factor +real xc, yc # computed centers +real xerr, yerr # estimate of centering error + +pointer sp, xm, ym + +begin + call smark (sp) + call salloc (xm, nx, TY_REAL) + call salloc (ym, ny, TY_REAL) + + # Compute the marginal distributions. + call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny) + call adivkr (Memr[xm], real (ny), Memr[xm], nx) + call adivkr (Memr[ym], real (nx), Memr[ym], ny) + + # Get the centers and errors. + call ap_cmarg (Memr[xm], nx, norm, xc, xerr) + call ap_cmarg (Memr[ym], ny, norm, yc, yerr) + + call sfree (sp) + return (AP_OK) +end + + +# AP_CMARG -- Compute the center and estimate its error given the +# marginal distribution and the number of points. + +procedure ap_cmarg (a, npts, norm, xc, err) + +real a[npts] # array +int npts # number of points +real norm # the normalization factor +real xc # center value +real err # error + +int i, npos +real val, sumi, sumix, sumix2 + +begin + # Initialize. + npos = 0 + sumi = 0.0 + sumix = 0.0 + sumix2 = 0.0 + + # Accumulate the sums. + do i = 1, npts { + val = a[i] + if (val > 0.0) + npos = npos + 1 + sumi = sumi + val + sumix = sumix + val * i + sumix2 = sumix2 + val * i ** 2 + } + + # Compute the position and the error. + if (npos <= 0) { + xc = (1.0 + npts) / 2.0 + err = INDEFR + } else { + xc = sumix / sumi + err = (sumix2 / sumi - xc * xc) + if (err <= 0.0) { + err = 0.0 + } else { + err = sqrt (err / (sumi * norm)) + if (err > real (npts)) + err = INDEFR + } + } +end + + +# AP_MKMARG -- Accumulate the marginal distributions. + +procedure ap_mkmarg (ctrpix, xm, ym, nx, ny) + +real ctrpix[nx, ny] # pixels +real xm[nx] # x marginal distribution +real ym[ny] # y marginal distribution +int nx, ny # dimensions of array + +int i, j +real sum + +begin + # Compute the x marginal. + do i = 1, nx { + sum = 0.0 + do j = 1, ny + sum = sum + ctrpix[i,j] + xm[i] = sum + } + + # Compute the y marginal. + do j = 1, ny { + sum = 0.0 + do i = 1, nx + sum = sum + ctrpix[i,j] + ym[j] = sum + } +end diff --git a/noao/digiphot/apphot/center/apctrbuf.x b/noao/digiphot/apphot/center/apctrbuf.x new file mode 100644 index 00000000..0a8053b8 --- /dev/null +++ b/noao/digiphot/apphot/center/apctrbuf.x @@ -0,0 +1,115 @@ +include <imhdr.h> +include <math.h> +include <mach.h> +include "../lib/apphotdef.h" +include "../lib/centerdef.h" +include "../lib/center.h" + +# APCTRBUF -- Procedure to fetch the center pixels given the pointer to the +# IRAF image, the coordinates of the initial center and the width of the APPHOT +# centering box. + +int procedure apctrbuf (ap, im, wx, wy) + +pointer ap # pointer to apphot structure +pointer im # pointer to the IRAF image +real wx, wy # center coordinates + +int icpix +pointer ctr +real cpix, gdatamin, gdatamax, datamin, datamax +pointer ap_ctrpix() + +begin + # Get pointer to centering structure. + ctr = AP_PCENTER(ap) + + # Check for 0 sized aperture. + if (AP_CAPERT(ctr) <= 0.0) + return (AP_CTR_NOAREA) + + # Get the centering buffer of pixels. + cpix = max (1.0, AP_CAPERT(ctr) * AP_SCALE(ap)) + icpix = 2 * int (cpix) + 1 + if (AP_CTRPIX(ctr) != NULL) + call mfree (AP_CTRPIX(ctr), TY_REAL) + AP_CTRPIX(ctr) = ap_ctrpix (im, wx, wy, icpix, + AP_CXC(ctr), AP_CYC(ctr), AP_CNX(ctr), AP_CNY(ctr)) + if (AP_CTRPIX(ctr) == NULL) + return (AP_CTR_NOAREA) + + # Compute the data limits. + if (IS_INDEFR(AP_DATAMIN(ap))) + gdatamin = -MAX_REAL + else + gdatamin = AP_DATAMIN(ap) + if (IS_INDEFR(AP_DATAMAX(ap))) + gdatamax = MAX_REAL + else + gdatamax = AP_DATAMAX(ap) + call alimr (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) * AP_CNY(ctr), + datamin, datamax) + + if (datamin < gdatamin || datamax > gdatamax) + return (AP_CTR_BADDATA) + else if (AP_CNX(ctr) < icpix || AP_CNY(ctr) < icpix) + return (AP_CTR_OUTOFBOUNDS) + else + return (AP_OK) +end + + +# AP_CTRPIX -- Procedure to fetch the pixels to be used for centering. + +pointer procedure ap_ctrpix (im, wx, wy, capert, xc, yc, nx, ny) + +pointer im # pointer to IRAF image +real wx, wy # center of subraster to be extracted +int capert # width of subraster to be extracted +real xc, yc # center of extracted subraster +int nx, ny # dimensions of extracted subraster + +int i, ncols, nlines, c1, c2, l1, l2, half_capert +pointer buf, lbuf +real xc1, xc2, xl1, xl2 +pointer imgs2r() + +begin + # Check for nonsensical input. + half_capert = (capert - 1) / 2 + if (half_capert <= 0) + return (NULL) + + # Test for out of bounds pixels + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + xc1 = wx - half_capert + xc2 = wx + half_capert + xl1 = wy - half_capert + xl2 = wy + half_capert + if (xc1 > real (ncols) || xc2 < 1.0 || xl1 > real (nlines) || xl2 < 1.0) + return (NULL) + + # Get column and line limits, dimensions and center of subraster. + c1 = max (1.0, min (real (ncols), xc1)) + 0.5 + c2 = min (real (ncols), max (1.0, xc2)) + 0.5 + l1 = max (1.0, min (real (nlines), xl1)) + 0.5 + l2 = min (real (nlines), max (1.0, xl2)) + 0.5 + nx = c2 - c1 + 1 + ny = l2 - l1 + 1 + xc = wx - c1 + 1 + yc = wy - l1 + 1 + + # Get pixels and return. + if (nx < 1 && ny < 1) + return (NULL) + else { + call malloc (buf, nx * ny, TY_REAL) + lbuf = buf + do i = l1, l2 { + call amovr (Memr[imgs2r (im, c1, c2, i, i)], Memr[lbuf], nx) + lbuf = lbuf + nx + } + return (buf) + } +end diff --git a/noao/digiphot/apphot/center/apfitcen.x b/noao/digiphot/apphot/center/apfitcen.x new file mode 100644 index 00000000..bd7224e8 --- /dev/null +++ b/noao/digiphot/apphot/center/apfitcen.x @@ -0,0 +1,291 @@ +include <mach.h> +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/noisedef.h" +include "../lib/centerdef.h" +include "../lib/center.h" + +define CONVERT .424660900 # conversion from fwhmpsf to sigma + +# APFITCENTER -- Procedure to fit the centers using either 1) the intensity +# weighted mean of the marginal distributions, 2) a Gaussian fit to the +# marginals assuming a fixed sigma or 3) a simplified version of the optimal +# filtering method using a Gaussian of fixed sigma to model the profile. + +int procedure apfitcenter (ap, im, wx, wy) + +pointer ap # pointer to the apphot structure +pointer im # pointer to the IRAF image +real wx, wy # object coordinates + +int i, niter, ier, fier, lowsnr +pointer ctr, nse +real cthreshold, datamin, datamax, owx, owy, xshift, yshift, med +int apctrbuf(), ap_ctr1d(), ap_mctr1d(), ap_gctr1d(), ap_lgctr1d() +real ap_csnratio(), amedr() + +begin + ctr = AP_PCENTER(ap) + nse = AP_NOISE(ap) + ier = AP_OK + + # Initialize. + AP_CXCUR(ctr) = wx + AP_CYCUR(ctr) = wy + AP_OXINIT(ctr) = INDEFR + AP_OYINIT(ctr) = INDEFR + AP_XCENTER(ctr) = INDEFR + AP_YCENTER(ctr) = INDEFR + AP_OXCENTER(ctr) = INDEFR + AP_OYCENTER(ctr) = INDEFR + AP_XSHIFT(ctr) = 0.0 + AP_YSHIFT(ctr) = 0.0 + AP_OXSHIFT(ctr) = 0.0 + AP_OYSHIFT(ctr) = 0.0 + AP_XERR(ctr) = INDEFR + AP_YERR(ctr) = INDEFR + + # Return input coordinates if centering is disabled. + if (IS_INDEFR(wx) || IS_INDEFR(wy)) { + return (AP_CTR_NOAREA) + } else if (AP_CENTERFUNCTION(ctr) == AP_NONE) { + AP_XCENTER(ctr) = wx + AP_YCENTER(ctr) = wy + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1) + call ap_ltoo (ap, wx, wy, AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1) + case WCS_TV: + call ap_ltov (im, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1) + call ap_ltov (im, wx, wy, AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1) + default: + AP_OXINIT(ctr) = wx + AP_OYINIT(ctr) = wy + AP_OXCENTER(ctr) = wx + AP_OYCENTER(ctr) = wy + } + return (AP_OK) + } + + # Intialize. + owx = wx + owy = wy + niter = 0 + if (IS_INDEFR(AP_SKYSIGMA(nse)) || IS_INDEFR(AP_CTHRESHOLD(ctr)) || + AP_CTHRESHOLD(ctr) <= 0.0) + cthreshold = 0.0 + else + cthreshold = AP_CTHRESHOLD(ctr) * AP_SKYSIGMA(nse) + + repeat { + + # Set initial cursor position. + AP_CXCUR(ctr) = owx + AP_CYCUR(ctr) = owy + + # Get centering pixels. + ier = apctrbuf (ap, im, owx, owy) + if (ier == AP_CTR_NOAREA) { + AP_XCENTER(ctr) = wx + AP_YCENTER(ctr) = wy + AP_XSHIFT(ctr) = 0.0 + AP_YSHIFT(ctr) = 0.0 + AP_XERR(ctr) = INDEFR + AP_YERR(ctr) = INDEFR + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1) + call ap_ltoo (ap, wx, wy, AP_OXCENTER(ctr), + AP_OYCENTER(ctr), 1) + case WCS_TV: + call ap_ltov (im, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1) + call ap_ltov (im, wx, wy, AP_OXCENTER(ctr), + AP_OYCENTER(ctr), 1) + default: + AP_OXINIT(ctr) = wx + AP_OYINIT(ctr) = wy + AP_OXCENTER(ctr) = wx + AP_OYCENTER(ctr) = wy + } + AP_OXSHIFT(ctr) = 0.0 + AP_OYSHIFT(ctr) = 0.0 + return (ier) + } + + # Clean the subraster. + if (AP_CLEAN(ctr) == YES) + call apclean (ap, Memr[AP_CTRPIX(ctr)], AP_CNX(ctr), + AP_CNY(ctr), AP_CXC(ctr), AP_CYC(ctr)) + + # Compute the datalimits. + if (cthreshold <= 0.0) + call alimr (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) * AP_CNY(ctr), + datamin, datamax) + else { + datamin = MAX_REAL + datamax = -MAX_REAL + do i = 1, AP_CNY(ctr) { + med = amedr (Memr[AP_CTRPIX(ctr)+(i-1)*AP_CNX(ctr)], + AP_CNX(ctr)) + if (med < datamin) + datamin = med + if (med > datamax) + datamax = med + } + } + + # Apply threshold and check for positive or negative features. + if (AP_POSITIVE(ap) == YES) { + call apsetr (ap, CDATALIMIT, datamin) + call asubkr (Memr[AP_CTRPIX(ctr)], datamin + + cthreshold, Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) * + AP_CNY(ctr)) + call amaxkr (Memr[AP_CTRPIX(ctr)], 0.0, Memr[AP_CTRPIX(ctr)], + AP_CNX(ctr) * AP_CNY(ctr)) + } else { + call apsetr (ap, CDATALIMIT, datamax) + call anegr (Memr[AP_CTRPIX(ctr)], Memr[AP_CTRPIX(ctr)], + AP_CNX(ctr) * AP_CNY(ctr)) + call aaddkr (Memr[AP_CTRPIX(ctr)], datamax - + cthreshold, Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) * + AP_CNY(ctr)) + call amaxkr (Memr[AP_CTRPIX(ctr)], 0.0, + Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) * AP_CNY(ctr)) + } + + # Test signal to noise ratio. + if (ap_csnratio (ap, Memr[AP_CTRPIX(ctr)], AP_CNX(ctr), + AP_CNY(ctr), 0.0) < AP_MINSNRATIO(ctr)) + lowsnr = YES + else + lowsnr = NO + + # Compute the x and y centers. + switch (AP_CENTERFUNCTION(ctr)) { + + case AP_CENTROID1D: + if (AP_CTHRESHOLD(ctr) <= 0.0) { + fier = ap_mctr1d (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr), + AP_CNY(ctr), AP_EPADU(nse), AP_XCENTER(ctr), + AP_YCENTER(ctr), AP_XERR(ctr), AP_YERR(ctr)) + if (IS_INDEFR (AP_XERR(ctr))) + AP_XCENTER(ctr) = AP_CXC(ctr) + if (IS_INDEFR (AP_YERR(ctr))) + AP_YCENTER(ctr) = AP_CYC(ctr) + } else { + fier = ap_ctr1d (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr), + AP_CNY(ctr), AP_EPADU(nse), AP_XCENTER(ctr), + AP_YCENTER(ctr), AP_XERR(ctr), AP_YERR(ctr)) + if (IS_INDEFR (AP_XERR(ctr))) + AP_XCENTER(ctr) = AP_CXC(ctr) + if (IS_INDEFR (AP_YERR(ctr))) + AP_YCENTER(ctr) = AP_CYC(ctr) + } + + case AP_GAUSS1D: + + fier = ap_gctr1d (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr), + AP_CNY(ctr), CONVERT * AP_FWHMPSF(ap) * AP_SCALE(ap), + AP_CMAXITER(ctr), AP_XCENTER(ctr), AP_YCENTER(ctr), + AP_XERR(ctr), AP_YERR(ctr)) + + case AP_OFILT1D: + + fier = ap_lgctr1d (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr), + AP_CNY(ctr), AP_CXC(ctr), AP_CYC(ctr), CONVERT * + AP_FWHMPSF(ap) * AP_SCALE(ap), AP_CMAXITER(ctr), + AP_EPADU(nse), AP_SKYSIGMA(nse), AP_XCENTER(ctr), + AP_YCENTER(ctr), AP_XERR(ctr), AP_YERR(ctr)) + + default: + # do nothing gracefully + } + + # Confine the next x center to the data box. + AP_XCENTER(ctr) = max (0.5, min (AP_CNX(ctr) + 0.5, + AP_XCENTER(ctr))) + xshift = AP_XCENTER(ctr) - AP_CXC(ctr) + AP_XCENTER(ctr) = xshift + owx + AP_XSHIFT(ctr) = AP_XCENTER(ctr) - wx + + # Confine the next y center to the data box. + AP_YCENTER(ctr) = max (0.5, min (AP_CNY(ctr) + 0.5, + AP_YCENTER(ctr))) + yshift = AP_YCENTER(ctr) - AP_CYC(ctr) + AP_YCENTER(ctr) = yshift + owy + AP_YSHIFT(ctr) = AP_YCENTER(ctr) - wy + + switch (AP_WCSOUT(ap)) { + case WCS_PHYSICAL, WCS_WORLD: + call ap_ltoo (ap, AP_XCENTER(ctr), AP_YCENTER(ctr), + AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1) + call ap_ltoo (ap, AP_XCENTER(ctr) - AP_XSHIFT(ctr), + AP_YCENTER(ctr) - AP_YSHIFT(ctr), AP_OXINIT(ctr), + AP_OYINIT(ctr), 1) + AP_OXSHIFT(ctr) = AP_OXCENTER(ctr) - AP_OXINIT(ctr) + AP_OYSHIFT(ctr) = AP_OYCENTER(ctr) - AP_OYINIT(ctr) + case WCS_TV: + call ap_ltov (im, AP_XCENTER(ctr), AP_YCENTER(ctr), + AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1) + call ap_ltov (im, AP_XCENTER(ctr) - AP_XSHIFT(ctr), + AP_YCENTER(ctr) - AP_YSHIFT(ctr), AP_OXINIT(ctr), + AP_OYINIT(ctr), 1) + AP_OXSHIFT(ctr) = AP_OXCENTER(ctr) - AP_OXINIT(ctr) + AP_OYSHIFT(ctr) = AP_OYCENTER(ctr) - AP_OYINIT(ctr) + default: + AP_OXINIT(ctr) = AP_XCENTER(ctr) - AP_XSHIFT(ctr) + AP_OYINIT(ctr) = AP_YCENTER(ctr) - AP_YSHIFT(ctr) + AP_OXCENTER(ctr) = AP_XCENTER(ctr) + AP_OYCENTER(ctr) = AP_YCENTER(ctr) + AP_OXSHIFT(ctr) = AP_XSHIFT(ctr) + AP_OYSHIFT(ctr) = AP_YSHIFT(ctr) + + } + + # Setup for next iteration. + niter = niter + 1 + owx = AP_XCENTER(ctr) + owy = AP_YCENTER(ctr) + + } until ((fier != AP_OK && fier != AP_CTR_NOCONVERGE) || + (niter >= AP_CMAXITER(ctr)) || (abs (xshift) < 1.0 && + abs (yshift) < 1.0)) + + # Return appropriate error code. + if (fier != AP_OK) { + AP_XCENTER(ctr) = wx + AP_YCENTER(ctr) = wy + AP_XSHIFT(ctr) = 0.0 + AP_YSHIFT(ctr) = 0.0 + AP_XERR(ctr) = INDEFR + AP_YERR(ctr) = INDEFR + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1) + call ap_ltoo (ap, wx, wy, AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1) + case WCS_TV: + call ap_ltov (im, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1) + call ap_ltov (im, wx, wy, AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1) + default: + AP_OXINIT(ctr) = wx + AP_OYINIT(ctr) = wy + AP_OXCENTER(ctr) = wx + AP_OYCENTER(ctr) = wy + } + AP_OXSHIFT(ctr) = 0.0 + AP_OYSHIFT(ctr) = 0.0 + return (fier) + } else if (ier == AP_CTR_BADDATA) { + return (AP_CTR_BADDATA) + } else if (lowsnr == YES) { + return (AP_CTR_LOWSNRATIO) + } else if (abs (AP_XSHIFT(ctr)) > (AP_MAXSHIFT(ctr) * AP_SCALE(ap))) { + return (AP_CTR_BADSHIFT) + } else if (abs (AP_YSHIFT(ctr)) > (AP_MAXSHIFT(ctr) * AP_SCALE(ap))) { + return (AP_CTR_BADSHIFT) + } else if (ier == AP_CTR_OUTOFBOUNDS) { + return (AP_CTR_OUTOFBOUNDS) + } else { + return (AP_OK) + } +end diff --git a/noao/digiphot/apphot/center/apgcpars.x b/noao/digiphot/apphot/center/apgcpars.x new file mode 100644 index 00000000..e40db12e --- /dev/null +++ b/noao/digiphot/apphot/center/apgcpars.x @@ -0,0 +1,27 @@ +include "../lib/display.h" +include "../lib/noise.h" +include "../lib/center.h" + +# AP_GCPARS -- Procedure to read in the centering parameters from the +# appropriate parameters files. + +procedure ap_gcpars (ap) + +pointer ap # pointer to the apphot structure + +bool clgetb() +int btoi() + +begin + # Initialize the structure. + call apcinit (ap, AP_CENTROID1D, 2.5, 2.0, AP_NPOISSON) + + # Get the data dependent parameters. + call ap_gdapars (ap) + + # Get the centering algorithm parameters. + call ap_gcepars (ap) + + # Make radial plots on stdgraph. + call apseti (ap, RADPLOTS, btoi (clgetb ("radplots"))) +end diff --git a/noao/digiphot/apphot/center/apgctr1d.x b/noao/digiphot/apphot/center/apgctr1d.x new file mode 100644 index 00000000..e0ac8077 --- /dev/null +++ b/noao/digiphot/apphot/center/apgctr1d.x @@ -0,0 +1,105 @@ +include <math/nlfit.h> +include "../lib/center.h" + +define NPARS 4 # the total number of parameters +define NAPARS 3 # the total number of active parameters +define TOL 0.001 # the tolerance for convergence + +# AP_GCTR1D -- Procedure to compute the x and y centers from the 1D marginal +# distributions using 1D Gaussian fits. Three parameters are fit for each +# marginal, the amplitude, the center of the Gaussian function itself +# and a constant background value. The sigma is set by the user and is +# assumed to be fixed. + +int procedure ap_gctr1d (ctrpix, nx, ny, sigma, maxiter, xc, yc, xerr, yerr) + +real ctrpix[nx, ny] # data subarray +int nx, ny # dimensions of data subarray +real sigma # sigma of PSF +int maxiter # maximum number of iterations +real xc, yc # computed centers +real xerr, yerr # estimate of centering error + +extern cgauss1d, cdgauss1d +int i, minel, maxel, xier, yier, npar, npts +pointer sp, x, xm, ym, w, fit, list, nl +real chisqr, variance, p[NPARS], dp[NPARS] +int locpr() + +begin + # Check the number of points. + if (nx < NAPARS || ny < NAPARS) + return (AP_CTR_NTOO_SMALL) + npts = max (nx, ny) + + call smark (sp) + call salloc (list, NAPARS, TY_INT) + call salloc (xm, nx, TY_REAL) + call salloc (ym, ny, TY_REAL) + call salloc (x, npts, TY_REAL) + call salloc (w, npts, TY_REAL) + call salloc (fit, npts, TY_REAL) + do i = 1, npts + Memr[x+i-1] = i + + # Compute the marginal distributions. + call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny) + call adivkr (Memr[xm], real (nx), Memr[xm], nx) + call adivkr (Memr[ym], real (ny), Memr[ym], ny) + + # Specify which parameters are to be fit. + Memi[list] = 1 + Memi[list+1] = 2 + Memi[list+2] = 4 + + # Initialize the x fit parameters. + call ap_alimr (Memr[xm], nx, p[4], p[1], minel, maxel) + p[1] = p[1] - p[4] + p[2] = maxel + p[3] = sigma ** 2 + + # Compute the x center and error. + call nlinitr (nl, locpr (cgauss1d), locpr (cdgauss1d), p, dp, NPARS, + Memi[list], NAPARS, TOL, maxiter) + call nlfitr (nl, Memr[x], Memr[xm], Memr[w], nx, 1, WTS_UNIFORM, xier) + call nlvectorr (nl, Memr[x], Memr[fit], nx, 1) + call nlpgetr (nl, p, npar) + call nlerrorsr (nl, Memr[xm], Memr[fit], Memr[w], nx, variance, + chisqr, dp) + call nlfreer (nl) + xc = p[2] + xerr = dp[2] / sqrt (real (nx)) + if (xerr > real (nx)) + xerr = INDEFR + + # Initialize the y fit parameters. + call ap_alimr (Memr[ym], ny, p[4], p[1], minel, maxel) + p[1] = p[1] - p[4] + p[2] = maxel + p[3] = sigma ** 2 + + # Fit the y marginal. + call nlinitr (nl, locpr (cgauss1d), locpr (cdgauss1d), p, dp, NPARS, + Memi[list], NAPARS, TOL, maxiter) + call nlfitr (nl, Memr[x], Memr[ym], Memr[w], ny, 1, WTS_UNIFORM, yier) + call nlvectorr (nl, Memr[x], Memr[fit], ny, 1) + call nlpgetr (nl, p, npar) + call nlerrorsr (nl, Memr[ym], Memr[fit], Memr[w], ny, variance, + chisqr, dp) + call nlfreer (nl) + yc = p[2] + yerr = dp[2] / sqrt (real (ny)) + if (yerr > real (ny)) + yerr = INDEFR + + # Return the appropriate error code. + call sfree (sp) + if (xier == NO_DEG_FREEDOM || yier == NO_DEG_FREEDOM) + return (AP_CTR_NTOO_SMALL) + else if (xier == SINGULAR || yier == SINGULAR) + return (AP_CTR_SINGULAR) + else if (xier == NOT_DONE || yier == NOT_DONE) + return (AP_CTR_NOCONVERGE) + else + return (AP_OK) +end diff --git a/noao/digiphot/apphot/center/apictr.x b/noao/digiphot/apphot/center/apictr.x new file mode 100644 index 00000000..3a3bd31e --- /dev/null +++ b/noao/digiphot/apphot/center/apictr.x @@ -0,0 +1,48 @@ +# AP_ICTR -- Given a subraster of pixels and an initial center compute +# the centroid of the pixels. + +procedure ap_ictr (im, wx, wy, radius, emission, xcenter, ycenter) + +pointer im # pointer to the iraf image +real wx # initial x coordinate +real wy # initial y coordinate +int radius # half width of centering box +int emission # emission feature +real xcenter # fitted x coordinate +real ycenter # fitted y coordinate + +int nx, ny, ier +pointer cbuf +real xc, yc, datamin, datamax, junk +int ap_mctr1d() +pointer ap_ctrpix() + +begin + # Get the pixels. + cbuf = ap_ctrpix (im, wx, wy, radius, xc, yc, nx, ny) + if (cbuf == NULL) { + xcenter = wx + ycenter = wy + return + } + + # Fit the center of the subraster. + call alimr (Memr[cbuf], nx * ny, datamin, datamax) + if (emission == YES) + call asubkr (Memr[cbuf], datamin, Memr[cbuf], nx * ny) + else { + call anegr (Memr[cbuf], Memr[cbuf], nx * ny) + call aaddkr (Memr[cbuf], datamax, Memr[cbuf], nx * ny) + } + ier = ap_mctr1d (Memr[cbuf], nx, ny, 1.0, xcenter, ycenter, junk, junk) + + # Compute the center in image coordinates. + if (IS_INDEFR(xcenter)) + xcenter = wx + else + xcenter = xcenter + wx - xc + if (IS_INDEFR(ycenter)) + ycenter = wy + else + ycenter = ycenter + wy - yc +end diff --git a/noao/digiphot/apphot/center/aplgctr1d.x b/noao/digiphot/apphot/center/aplgctr1d.x new file mode 100644 index 00000000..ce6d7c53 --- /dev/null +++ b/noao/digiphot/apphot/center/aplgctr1d.x @@ -0,0 +1,87 @@ +include <math.h> +include "../lib/center.h" + +define TOL 0.001 # tolerance for fitting algorithm + +# AP_LGCTR1D -- Procedure to compute the center from the 1D marginal +# distributions using a simplified version of the optimal filtering +# technique and addopting a Gaussian model for the fit. The method +# is streamlined by replacing the Gaussian with a simple triangle +# following L.Goad. + +int procedure ap_lgctr1d (ctrpix, nx, ny, cx, cy, sigma, maxiter, norm, + skysigma, xc, yc, xerr, yerr) + +real ctrpix[nx, ny] # object to be centered +int nx, ny # dimensions of subarray +real cx, cy # center in subraster coordinates +real sigma # sigma of PSF +int maxiter # maximum number of iterations +real norm # the normalization factor +real skysigma # standard deviation of the pixels +real xc, yc # computed centers +real xerr, yerr # first guess at errors + +int nxiter, nyiter +pointer sp, xm, ym +real ratio, constant + +int aptopt() +real asumr() + +begin + # Allocate working space. + call smark (sp) + call salloc (xm, nx, TY_REAL) + call salloc (ym, ny, TY_REAL) + + # Compute the marginal distributions. + call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny) + xerr = asumr (Memr[xm], nx) + yerr = asumr (Memr[ym], ny) + call adivkr (Memr[xm], real (nx), Memr[xm], nx) + call adivkr (Memr[ym], real (ny), Memr[ym], ny) + + # Compute the x center and error. + xc = cx + nxiter = aptopt (Memr[xm], nx, xc, sigma, TOL, maxiter, YES) + if (xerr <= 0.0) + xerr = INDEFR + else { + if (IS_INDEFR(skysigma)) + constant = 0.0 + else + constant = 4.0 * SQRTOFPI * sigma * skysigma ** 2 + ratio = constant / xerr + xerr = sigma ** 2 / (xerr * norm) + xerr = sqrt (max (xerr, ratio * xerr)) + if (xerr > real (nx)) + xerr = INDEFR + } + + # Compute the y center and error. + yc = cy + nyiter = aptopt (Memr[ym], ny, yc, sigma, TOL, maxiter, YES) + if (yerr <= 0.0) + yerr = INDEFR + else { + if (IS_INDEFR(skysigma)) + constant = 0.0 + else + constant = 4.0 * SQRTOFPI * sigma * skysigma ** 2 + ratio = constant / yerr + yerr = sigma ** 2 / (yerr * norm) + yerr = sqrt (max (yerr, ratio * yerr)) + if (yerr > real (ny)) + yerr = INDEFR + } + + # Return appropriate error code. + call sfree (sp) + if (nxiter < 0 || nyiter < 0) + return (AP_CTR_SINGULAR) + else if (nxiter > maxiter || nyiter > maxiter) + return (AP_CTR_NOCONVERGE) + else + return (AP_OK) +end diff --git a/noao/digiphot/apphot/center/apmctr1d.x b/noao/digiphot/apphot/center/apmctr1d.x new file mode 100644 index 00000000..5dc47f60 --- /dev/null +++ b/noao/digiphot/apphot/center/apmctr1d.x @@ -0,0 +1,88 @@ +include <mach.h> +include "../lib/center.h" + +# AP_MCTR1D -- Procedure to compute the center from the 1D marginal +# distributions. + +int procedure ap_mctr1d (ctrpix, nx, ny, norm, xc, yc, xerr, yerr) + +real ctrpix[nx, ny] # object to be centered +int nx, ny # dimensions of subarray +real norm # the normalization factor +real xc, yc # computed centers +real xerr, yerr # estimate of centering error + +pointer sp, xm, ym + +begin + call smark (sp) + call salloc (xm, nx, TY_REAL) + call salloc (ym, ny, TY_REAL) + + # Compute the marginal distributions. + call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny) + call adivkr (Memr[xm], real (ny), Memr[xm], nx) + call adivkr (Memr[ym], real (nx), Memr[ym], ny) + + # Get the centers and errors. + call ap_cmmarg (Memr[xm], nx, norm, xc, xerr) + call ap_cmmarg (Memr[ym], ny, norm, yc, yerr) + + call sfree (sp) + return (AP_OK) +end + + +# AP_CMMARG -- Compute the center and estimate its error given the +# marginal distribution and the number of points. The center is estimated +# using a centroid which is an optimal solution for a gaussian star with +# no background noise. The centering error is computing using a variance +# estimate based on a gaussian star with no background noise. + +procedure ap_cmmarg (a, npts, norm, xc, err) + +real a[npts] # array +int npts # number of points +real norm # the normalization factor +real xc # center value +real err # error + +int i, npos +real sumi, sumix, sumix2, mean, val +real asumr() + +begin + # Initialize. + mean = asumr (a, npts) / npts + npos = 0 + sumi = 0.0 + sumix = 0.0 + sumix2 = 0.0 + + # Accumulate the sums. + do i = 1, npts { + val = (a[i] - mean) + if (val > 0.0) { + npos = npos + 1 + sumi = sumi + val + sumix = sumix + val * i + sumix2 = sumix2 + val * i ** 2 + } + } + + # Compute the position and the error. + if (npos <= 0) { + xc = (1.0 + npts) / 2.0 + err = INDEFR + } else { + xc = sumix / sumi + err = (sumix2 / sumi - xc * xc) + if (err <= 0.0) { + err = 0.0 + } else { + err = sqrt (err / ((sumi + mean * npts) * norm)) + if (err > real (npts)) + err = INDEFR + } + } +end diff --git a/noao/digiphot/apphot/center/appcenter.x b/noao/digiphot/apphot/center/appcenter.x new file mode 100644 index 00000000..c2590fa3 --- /dev/null +++ b/noao/digiphot/apphot/center/appcenter.x @@ -0,0 +1,84 @@ +include "../lib/apphot.h" +include "../lib/center.h" + +# AP_PCENTER -- Procedure to write the output of centering task to a file. + +procedure ap_pcenter (ap, fd, id, lid, ier) + +pointer ap # pointer to apphot structure +int fd # output file descriptor +int id # id of the star +int lid # list number +int ier # comment string + +real xpos, ypos +real apstatr() + +begin + if (fd == NULL) + return + + # Print description of object. + xpos = apstatr (ap, OXINIT) + ypos = apstatr (ap, OYINIT) + call ap_wid (ap, fd, xpos, ypos, id, lid, '\\') + + # Print the computed centers. + call ap_wcres (ap, fd, ier, ' ') +end + + +# AP_QCENTER -- Procedure to print a quick summary of the center task output +# on the standard output. + +procedure ap_qcenter (ap, ier) + +pointer ap # pointer to apphot structure +int ier # comment string + +real owx, owy, wx, wy +pointer sp, imname +real apstatr() + +begin + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + + owx = apstatr (ap, OXINIT) + owy = apstatr (ap, OYINIT) + wx = apstatr (ap, OXCENTER) + wy = apstatr (ap, OYCENTER) + + call apstats (ap, IMROOT, Memc[imname], SZ_FNAME) + call printf ( "%s %7.2f %7.2f %7.2f %7.2f ") + call pargstr (Memc[imname]) + call pargr (owx) + call pargr (owy) + call pargr (wx) + call pargr (wy) + call printf ("%6.2f %6.2f %s\n") + call pargr (apstatr (ap, XERR)) + call pargr (apstatr (ap, YERR)) + if (ier != AP_OK) + call pargstr ("err") + else + call pargstr ("ok") + + call sfree (sp) +end + + +# AP_CPHDR -- Procedure to write the banner for the center task to the +# output file. + +procedure ap_cphdr (ap, fd) + +pointer ap # apphot descriptor +int fd # output file descriptor + +begin + if (fd == NULL) + return + call ap_idhdr (ap, fd) + call ap_chdr (ap, fd) +end diff --git a/noao/digiphot/apphot/center/appcpars.x b/noao/digiphot/apphot/center/appcpars.x new file mode 100644 index 00000000..272fcc07 --- /dev/null +++ b/noao/digiphot/apphot/center/appcpars.x @@ -0,0 +1,21 @@ +include "../lib/display.h" + +# AP_PCPARS -- Procedure to write out the current centering parameters +# to the current parameter files. + +procedure ap_pcpars (ap) + +pointer ap # pointer to apphot structure +bool itob() +int apstati() + +begin + # Write the data dependent parameters. + call ap_dapars (ap) + + # Write the centering parameters. + call ap_cepars (ap) + + # Set the plotting command. + call clputb ("radplots", itob (apstati (ap, RADPLOTS))) +end diff --git a/noao/digiphot/apphot/center/aprefitcen.x b/noao/digiphot/apphot/center/aprefitcen.x new file mode 100644 index 00000000..c50e2efb --- /dev/null +++ b/noao/digiphot/apphot/center/aprefitcen.x @@ -0,0 +1,193 @@ +include <mach.h> +include "../lib/apphotdef.h" +include "../lib/apphot.h" +include "../lib/noisedef.h" +include "../lib/centerdef.h" +include "../lib/center.h" + +define CONVERT .424660900 # conversion factor from fwhmpsf to sigma + +# APREFITCENTER -- Procedure to refit the centers assuming that the appropriate +# pixel buffer is in memory. See apfitcenter for further information. + +int procedure aprefitcenter (ap, im, ier) + +pointer ap # pointer to the apphot structure +pointer im # the input image descriptor +int ier # previous error code + +int fier +pointer cen, nse +int ap_ctr1d(), ap_mctr1d(), ap_gctr1d(), ap_lgctr1d() + +begin + cen = AP_PCENTER(ap) + nse = AP_NOISE(ap) + + # Initialize + AP_XCENTER(cen) = AP_CXCUR(cen) + AP_YCENTER(cen) = AP_CYCUR(cen) + AP_XSHIFT(cen) = 0.0 + AP_YSHIFT(cen) = 0.0 + AP_OXSHIFT(cen) = 0.0 + AP_OYSHIFT(cen) = 0.0 + AP_XERR(cen) = INDEFR + AP_YERR(cen) = INDEFR + + # Return if the center is undefined. + if (IS_INDEFR(AP_CXCUR(cen)) || IS_INDEFR(AP_CYCUR(cen))) { + AP_OXINIT(cen) = INDEFR + AP_OYINIT(cen) = INDEFR + AP_OXCENTER(cen) = INDEFR + AP_OYCENTER(cen) = INDEFR + return (AP_CTR_NOAREA) + } + + # Convert the coordinates. + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXINIT(cen), + AP_OYINIT(cen), 1) + call ap_ltoo (ap, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXCENTER(cen), + AP_OYCENTER(cen), 1) + case WCS_TV: + call ap_ltov (im, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXINIT(cen), + AP_OYINIT(cen), 1) + call ap_ltov (im, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXCENTER(cen), + AP_OYCENTER(cen), 1) + default: + AP_OXINIT(cen) = AP_CXCUR(cen) + AP_OYINIT(cen) = AP_CYCUR(cen) + AP_OXCENTER(cen) = AP_CXCUR(cen) + AP_OYCENTER(cen) = AP_CYCUR(cen) + } + + # Return input coordinates if no center fitting. + if (AP_CENTERFUNCTION(cen) == AP_NONE) + return (AP_OK) + + # Choose the centering algorithm. + switch (AP_CENTERFUNCTION(cen)) { + + case AP_CENTROID1D: + + # Compute the x and y centroids. + if (AP_CTHRESHOLD(cen) <= 0.0) { + fier = ap_mctr1d (Memr[AP_CTRPIX(cen)], AP_CNX(cen), + AP_CNY(cen), AP_EPADU(nse), AP_XCENTER(cen), + AP_YCENTER(cen), AP_XERR(cen), AP_YERR(cen)) + if (IS_INDEFR(AP_XERR(cen))) + AP_XCENTER(cen) = AP_CXC(cen) + if (IS_INDEFR(AP_YERR(cen))) + AP_YCENTER(cen) = AP_CYC(cen) + } else { + fier = ap_ctr1d (Memr[AP_CTRPIX(cen)], AP_CNX(cen), + AP_CNY(cen), AP_EPADU(nse), AP_XCENTER(cen), + AP_YCENTER(cen), AP_XERR(cen), AP_YERR(cen)) + if (IS_INDEFR(AP_XERR(cen))) + AP_XCENTER(cen) = AP_CXC(cen) + if (IS_INDEFR(AP_YERR(cen))) + AP_YCENTER(cen) = AP_CYC(cen) + } + AP_XCENTER(cen) = AP_XCENTER(cen) + AP_CXCUR(cen) - AP_CXC(cen) + AP_YCENTER(cen) = AP_YCENTER(cen) + AP_CYCUR(cen) - AP_CYC(cen) + AP_XSHIFT(cen) = AP_XCENTER(cen) - AP_CXCUR(cen) + AP_YSHIFT(cen) = AP_YCENTER(cen) - AP_CYCUR(cen) + + case AP_GAUSS1D: + + # Compute the 1D Gaussian x and y centers. + fier = ap_gctr1d (Memr[AP_CTRPIX(cen)], AP_CNX(cen), AP_CNY(cen), + CONVERT * AP_FWHMPSF(ap) * AP_SCALE(ap), AP_CMAXITER(cen), + AP_XCENTER(cen), AP_YCENTER(cen), AP_XERR(cen), AP_YERR(cen)) + AP_XCENTER(cen) = AP_XCENTER(cen) + AP_CXCUR(cen) - AP_CXC(cen) + AP_YCENTER(cen) = AP_YCENTER(cen) + AP_CYCUR(cen) - AP_CYC(cen) + AP_XSHIFT(cen) = AP_XCENTER(cen) - AP_CXCUR(cen) + AP_YSHIFT(cen) = AP_YCENTER(cen) - AP_CYCUR(cen) + + case AP_OFILT1D: + + # Compute the Goad 1D x and y centers. + fier = ap_lgctr1d (Memr[AP_CTRPIX(cen)], AP_CNX(cen), AP_CNY(cen), + AP_CXC(cen), AP_CYC(cen), CONVERT * AP_FWHMPSF(ap) * + AP_SCALE(ap), AP_CMAXITER(cen), AP_EPADU(nse), + AP_SKYSIGMA(nse), AP_XCENTER(cen), AP_YCENTER(cen), + AP_XERR(cen), AP_YERR(cen)) + AP_XCENTER(cen) = AP_XCENTER(cen) + AP_CXCUR(cen) - AP_CXC(cen) + AP_YCENTER(cen) = AP_YCENTER(cen) + AP_CYCUR(cen) - AP_CYC(cen) + AP_XSHIFT(cen) = AP_XCENTER(cen) - AP_CXCUR(cen) + AP_YSHIFT(cen) = AP_YCENTER(cen) - AP_CYCUR(cen) + + default: + + # do nothing gracefully + } + + # Convert the coordinates. + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, AP_XCENTER(cen), AP_YCENTER(cen), + AP_OXCENTER(cen), AP_OYCENTER(cen), 1) + call ap_ltoo (ap, AP_XCENTER(cen) - AP_XSHIFT(cen), + AP_YCENTER(cen) - AP_YSHIFT(cen), AP_OXINIT(cen), + AP_OYINIT(cen), 1) + AP_OXSHIFT(cen) = AP_OXCENTER(cen) - AP_OXINIT(cen) + AP_OYSHIFT(cen) = AP_OYCENTER(cen) - AP_OYINIT(cen) + case WCS_TV: + call ap_ltov (im, AP_XCENTER(cen), AP_YCENTER(cen), + AP_OXCENTER(cen), AP_OYCENTER(cen), 1) + call ap_ltov (im, AP_XCENTER(cen) - AP_XSHIFT(cen), + AP_YCENTER(cen) - AP_YSHIFT(cen), AP_OXINIT(cen), + AP_OYINIT(cen), 1) + AP_OXSHIFT(cen) = AP_OXCENTER(cen) - AP_OXINIT(cen) + AP_OYSHIFT(cen) = AP_OYCENTER(cen) - AP_OYINIT(cen) + default: + AP_OXCENTER(cen) = AP_XCENTER(cen) + AP_OYCENTER(cen) = AP_YCENTER(cen) + AP_OXINIT(cen) = AP_XCENTER(cen) - AP_XSHIFT(cen) + AP_OYINIT(cen) = AP_YCENTER(cen) - AP_YSHIFT(cen) + AP_OXSHIFT(cen) = AP_XSHIFT(cen) + AP_OYSHIFT(cen) = AP_YSHIFT(cen) + } + + # Return appropriate error code. + if (fier != AP_OK) { + AP_XCENTER(cen) = AP_CXCUR(cen) + AP_YCENTER(cen) = AP_CYCUR(cen) + AP_XSHIFT(cen) = 0.0 + AP_YSHIFT(cen) = 0.0 + AP_XERR(cen) = INDEFR + AP_YERR(cen) = INDEFR + switch (AP_WCSOUT(ap)) { + case WCS_WORLD, WCS_PHYSICAL: + call ap_ltoo (ap, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXINIT(cen), + AP_OYINIT(cen), 1) + call ap_ltoo (ap, AP_CXCUR(cen), AP_CYCUR(cen), + AP_OXCENTER(cen), AP_OYCENTER(cen), 1) + case WCS_TV: + call ap_ltov (im, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXINIT(cen), + AP_OYINIT(cen), 1) + call ap_ltov (im, AP_CXCUR(cen), AP_CYCUR(cen), + AP_OXCENTER(cen), AP_OYCENTER(cen), 1) + default: + AP_OXCENTER(cen) = AP_CXCUR(cen) + AP_OYCENTER(cen) = AP_CYCUR(cen) + AP_OXINIT(cen) = AP_CXCUR(cen) + AP_OYINIT(cen) = AP_CYCUR(cen) + } + AP_OXSHIFT(cen) = 0.0 + AP_OYSHIFT(cen) = 0.0 + return (fier) + } else if (ier == AP_CTR_BADDATA) { + return (AP_CTR_BADDATA) + } else if (ier == AP_CTR_LOWSNRATIO) { + return (AP_CTR_LOWSNRATIO) + } else if (abs (AP_XSHIFT(cen)) > (AP_MAXSHIFT(cen) * AP_SCALE(ap))) { + return (AP_CTR_BADSHIFT) + } else if (abs (AP_YSHIFT(cen)) > (AP_MAXSHIFT(cen) * AP_SCALE(ap))) { + return (AP_CTR_BADSHIFT) + } else if (ier == AP_CTR_OUTOFBOUNDS) { + return (AP_CTR_OUTOFBOUNDS) + } else + return (AP_OK) +end diff --git a/noao/digiphot/apphot/center/center.key b/noao/digiphot/apphot/center/center.key new file mode 100644 index 00000000..f71cd8d3 --- /dev/null +++ b/noao/digiphot/apphot/center/center.key @@ -0,0 +1,78 @@ + Interactive Keystroke Commands + +? Print help +: Colon commands +v Verify the critical parameters +w Save the current parameters +d Plot radial profile of current star +i Interactively set parameters using current star +f Fit center of current star +spbar Fit center of current star, output results +m Move to next star in coordinate list +n Center next star in coordinate list, output results +l Center remaining stars in coordinate list, output results +e Print error messages +r Rewind the coordinate list +q Exit task + + + Colon Commands + +:show [data/center] List the parameters +:m [n] Move to next [nth] star in coordinate list +:n [n] Center next [nth] star in coordinate list, output results + + + Colon Parameter Editing Commands + +# Image and file name parameters + +:image [string] Image name +:coords [string] Coordinate file name +:output [string] Output file name + +# Data dependent parameters + +:scale [value] Image scale (units per pixel) +:fwhmpsf [value] Full-width half-maximum of PSF (scale units) +:emission [y/n] Emission feature (y), absorption (n) +:sigma [value] Standard deviation of sky (counts) +:datamin [value] Minimum good data value (counts) +:datamax [value] Maximum good data value (counts) + +# Noise parameters + +:noise [string] Noise model (constant|poisson) +:gain [string] Gain image header keyword +:ccdread [string] Readout noise image header keyword +:epadu [value] Gain (electrons per adu) +:readnoise [value] Readout noise (electrons) + +# Observations parameters + +:exposure [string] Exposure time image header keyword +:airmass [string] Airmass image header keyword +:filter [string] Filter image header keyword +:obstime [string] Time of observation image header keyword +:itime [value] Exposure time (time units) +:xairmass [value] Airmass value (number) +:ifilter [string] Filter id string +:otime [string] Time of observation (time units) + +# Centering parameters + +:calgorithm [string] Centering algorithm +:cbox [value] Width of centering box (scale units) +:cthreshold [value] Centering intensity threshold (sigma) +:cmaxiter [value] Maximum number of iterations +:maxshift [value] Maximum center shift (scale units) +:minsnratio [value] Minimum signal to noise for centering +:clean [y/n] Clean subraster before centering +:rclean [value] Cleaning radius (scale units) +:rclip [value] Clipping radius (scale units) +:kclean [value] Clean K-sigma rejection limit (sigma) + +# Plotting and marking parameters + +:mkcenter [y/n] Mark computed centers on the display +:radplot [y/n] Plot radial profile of object diff --git a/noao/digiphot/apphot/center/icenter.key b/noao/digiphot/apphot/center/icenter.key new file mode 100644 index 00000000..0e5379f0 --- /dev/null +++ b/noao/digiphot/apphot/center/icenter.key @@ -0,0 +1,12 @@ + Interactive Center Setup Menu + + v Mark and verify the critical center parameters (f,s,c) + + f Mark and verify the full-width half-maximum of the psf + s Mark and verify the standard deviation of the background + l Mark and verify the minimum good data value + u Mark and verify the maximum good data value + + c Mark and verify the centering box half-width + n Mark and verify the cleaning radius + p Mark and verify the clipping radius diff --git a/noao/digiphot/apphot/center/mkpkg b/noao/digiphot/apphot/center/mkpkg new file mode 100644 index 00000000..dde81929 --- /dev/null +++ b/noao/digiphot/apphot/center/mkpkg @@ -0,0 +1,54 @@ +# CENTER task + +$checkout libpkg.a ".." +$update libpkg.a +$checkin libpkg.a ".." +$exit + +libpkg.a: + apbcenter.x ../lib/apphot.h ../lib/display.h \ + <fset.h> + apcconfirm.x ../lib/apphot.h ../lib/noise.h \ + ../lib/center.h + apcencolon.x ../lib/apphot.h ../lib/center.h \ + ../lib/display.h ../lib/noise.h + apcenter.x ../lib/apphot.h ../lib/center.h \ + ../lib/display.h <ctype.h> \ + <gset.h> <imhdr.h> + apcerrors.x ../lib/center.h + apcfree.x ../lib/apphotdef.h ../lib/centerdef.h + apictr.x + apcinit.x ../lib/apphotdef.h ../lib/center.h \ + ../lib/centerdef.h + apclean.x ../lib/apphot.h ../lib/center.h \ + ../lib/noise.h + appcpars.x ../lib/display.h + apcplot.x ../lib/apphot.h ../lib/apphotdef.h \ + ../lib/center.h ../lib/centerdef.h \ + ../lib/noise.h <gset.h> \ + <pkg/gtools.h> <mach.h> + apcpshow.x ../lib/center.h ../lib/display.h + apcshow.x ../lib/display.h + apcsnratio.x ../lib/center.h ../lib/noise.h \ + <mach.h> + apctr1d.x ../lib/center.h <mach.h> + apctrbuf.x ../lib/apphotdef.h ../lib/center.h \ + ../lib/centerdef.h <imhdr.h> <math.h> \ + <mach.h> + apmctr1d.x ../lib/center.h <mach.h> + apfitcen.x ../lib/apphotdef.h ../lib/center.h \ + ../lib/centerdef.h ../lib/noisedef.h \ + ../lib/apphot.h <mach.h> + apgcpars.x ../lib/center.h ../lib/display.h \ + ../lib/noise.h + apgctr1d.x ../lib/center.h <math/nlfit.h> + aplgctr1d.x ../lib/center.h <math.h> + appcenter.x ../lib/apphot.h ../lib/center.h + apcradsetup.x ../lib/display.h + aprefitcen.x ../lib/apphotdef.h ../lib/center.h \ + ../lib/centerdef.h ../lib/noisedef.h \ + ../lib/apphot.h <mach.h> + t_center.x ../lib/apphot.h ../lib/noise.h \ + <fset.h> <gset.h> \ + <imhdr.h> + ; diff --git a/noao/digiphot/apphot/center/t_center.x b/noao/digiphot/apphot/center/t_center.x new file mode 100644 index 00000000..456b93d2 --- /dev/null +++ b/noao/digiphot/apphot/center/t_center.x @@ -0,0 +1,306 @@ +include <fset.h> +include <gset.h> +include <imhdr.h> +include "../lib/apphot.h" +include "../lib/noise.h" + +# T_CENTER -- Procedure to compute accurate centers for a list of objects +# given a list of starting centers, a set of starting parameters and a list +# of images. + +procedure t_center () + +pointer image # pointer to name of the image +pointer output # pointer to output file name +pointer coords # pointer to coordinate file +pointer plotfile # name of plot metacode file +pointer graphics # pointer to graphics device +pointer display # pointer to display device +int interactive # interactive mode +int cache # cache image buffer in memory +int verify # verify parameters +int update # update parameters +int verbose # verbose mode + +pointer sp, str, cname, outfname, ap, im, id, gd, mgd +int limlist, lclist, lolist, sid, lid, cl, pfd, out, root, stat, memstat +int imlist, clist, olist, wcs, req_size, old_size, buf_size + +pointer immap(), gopen() +int imtlen(), imtgetim(), clplen(), clgfil(), btoi(), fnldir(), strncmp() +int open(), strlen(), apcenter(), imtopenp(), clpopnu(), clgwrd() +int ap_memstat(), sizeof() +bool clgetb(), streq() +errchk gopen + +begin + # Allocate working space. + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (coords, SZ_FNAME, TY_CHAR) + call salloc (plotfile, SZ_FNAME, TY_CHAR) + call salloc (graphics, SZ_FNAME, TY_CHAR) + call salloc (display, SZ_FNAME, TY_CHAR) + call salloc (outfname, SZ_FNAME, TY_CHAR) + call salloc (cname, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Set standard output to flush on newline. + call fseti (STDOUT, F_FLUSHNL, YES) + + # Get input and output file names. + imlist = imtopenp ("image") + limlist = imtlen (imlist) + clist = clpopnu ("coords") + lclist = clplen (clist) + olist = clpopnu ("output") + lolist = clplen (olist) + + # Check that image and coordinate list lengths match. + if (limlist < 1 || (lclist > 1 && lclist != limlist)) { + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + call error (0, "Imcompatible image and coordinate list lengths") + } + + # Check that image and output list lengths match. + if (lolist > 1 && lolist != limlist) { + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + call error (0, "Imcompatible image and output list lengths") + } + + call clgstr ("icommands.p_filename", Memc[cname], SZ_FNAME) + if (Memc[cname] != EOS) + interactive = NO + #else if (lclist == 0) + #interactive = YES + else + interactive = btoi (clgetb ("interactive")) + cache = btoi (clgetb ("cache")) + verbose = btoi (clgetb ("verbose")) + verify = btoi (clgetb ("verify")) + update = btoi (clgetb ("update")) + + # Get the centering parameters. + call ap_gcpars (ap) + + # confirm the centering algorithm parameters. + if (verify == YES && interactive == NO) { + call ap_cconfirm (ap, NULL, 1) + if (update == YES) + call ap_pcpars (ap) + } + + # Get the wcs information. + wcs = clgwrd ("wcsin", Memc[str], SZ_LINE, WCSINSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the input coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call apseti (ap, WCSIN, wcs) + wcs = clgwrd ("wcsout", Memc[str], SZ_LINE, WCSOUTSTR) + if (wcs <= 0) { + call eprintf ( + "Warning: Setting the output coordinate system to logical\n") + wcs = WCS_LOGICAL + } + call apseti (ap, WCSOUT, wcs) + + # Get the graphics and display devices. + call clgstr ("graphics", Memc[graphics], SZ_FNAME) + call clgstr ("display", Memc[display], SZ_FNAME) + + # Open the graphics and image deisplay devices. + if (interactive == YES) { + if (Memc[graphics] == EOS) + gd = NULL + else { + iferr { + gd = gopen (Memc[graphics], APPEND+AW_DEFER, STDGRAPH) + } then { + call eprintf ( + "Warning: Error opening the graphics device.\n") + gd = NULL + } + } + if (Memc[display] == EOS) + id = NULL + else if (streq (Memc[graphics], Memc[display])) + id = gd + else { + iferr { + id = gopen (Memc[display], APPEND, STDIMAGE) + } then { + call eprintf ( + "Warning: Graphics overlay not available for display device.\n") + id = NULL + } + } + } else { + gd = NULL + id = NULL + } + + # Open the plot metacode file. + call clgstr ("plotfile", Memc[plotfile], SZ_FNAME) + if (Memc[plotfile] == EOS) + pfd = NULL + else + pfd = open (Memc[plotfile], APPEND, BINARY_FILE) + if (pfd != NULL) + mgd = gopen (Memc[graphics], NEW_FILE, pfd) + else + mgd = NULL + + # Begin looping over the image list. + sid = 1 + while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) { + + # Open image. + im = immap (Memc[image], READ_ONLY, 0) + call apimkeys (ap, im, Memc[image]) + + # Set the image display viewport. + if ((id != NULL) && (id != gd)) + call ap_gswv (id, Memc[image], im, 4) + + # Cache the input image pixels. + req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) * + sizeof (IM_PIXTYPE(im)) + memstat = ap_memstat (cache, req_size, old_size) + if (memstat == YES) + call ap_pcache (im, INDEFI, buf_size) + + # Open the coordinate file; where coords is assumed to be a simple + # text file in which the x and y positions are in columns 1 and 2 + # respectively and all remaining fields are ignored. + + if (lclist <= 0) { + cl = NULL + call strcpy ("", Memc[outfname], SZ_FNAME) + } else { + stat = clgfil (clist, Memc[coords], SZ_FNAME) + root = fnldir (Memc[coords], Memc[outfname], SZ_FNAME) + if (strncmp ("default", Memc[coords+root], 7) == 0 || root == + strlen (Memc[coords])) { + call ap_inname (Memc[image], Memc[outfname], "coo", + Memc[outfname], SZ_FNAME) + lclist = limlist + cl = open (Memc[outfname], READ_ONLY, TEXT_FILE) + } else if (stat != EOF) { + call strcpy (Memc[coords], Memc[outfname], SZ_FNAME) + cl = open (Memc[outfname], READ_ONLY, TEXT_FILE) + } else { + call apstats (ap, CLNAME, Memc[outfname], SZ_FNAME) + call seek (cl, BOF) + } + } + call apsets (ap, CLNAME, Memc[outfname]) + call apfroot (Memc[outfname], Memc[str], SZ_LINE) + call apsets (ap, CLROOT, Memc[str]) + + # Open output text file; if output is "default", dir$default or + # a directory specification then the extension "ctr" is added to + # the root image name and a suitable version number is appended to + # the output name. If the output string is null then no output + # file is created. + + if (lolist == 0) { + out = NULL + call strcpy ("", Memc[outfname], SZ_FNAME) + } else { + stat = clgfil (olist, Memc[output], SZ_FNAME) + root = fnldir (Memc[output], Memc[outfname], SZ_FNAME) + if (strncmp ("default", Memc[output+root], 7) == 0 || root == + strlen (Memc[output])) { + call apoutname (Memc[image], Memc[outfname], "ctr", + Memc[outfname], SZ_FNAME) + out = open (Memc[outfname], NEW_FILE, TEXT_FILE) + lolist = limlist + } else if (stat != EOF) { + call strcpy (Memc[output], Memc[outfname], SZ_FNAME) + out = open (Memc[outfname], NEW_FILE, TEXT_FILE) + } else + call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME) + } + call apsets (ap, OUTNAME, Memc[outfname]) + + # Fit the centers. + if (interactive == NO) { + if (Memc[cname] != EOS) + stat = apcenter (ap, im, cl, NULL, mgd, NULL, out, sid, + NO, cache) + else if (cl != NULL) { + lid = 1 + call apbcenter (ap, im, cl, out, sid, lid, mgd, id, + verbose) + stat = NO + } else + stat = NO + } else + stat = apcenter (ap, im, cl, gd, mgd, id, out, sid, YES, cache) + + # Close up image, coordinates, and output file. + call imunmap (im) + if (cl != NULL) { + if (lclist > 1) + call close (cl) + } + if (out != NULL && lolist != 1) { + call close (out) + if (sid <= 1) { + call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME) + call delete (Memc[outfname]) + } + sid = 1 + } + + # Uncache memory. + call fixmem (old_size) + + if (stat == YES) + break + + } + + # If only one coordinate file for a list of images close it. + if (cl != NULL && lclist == 1) + call close (cl) + + # If only one output file defined for a list of images close it. + if (out != NULL && lolist == 1) { + call close (out) + if (sid <= 1) { + call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME) + call delete (Memc[outfname]) + } + } + + # Close up graphics and image display streams and the plot files. + if (id == gd && id != NULL) + call gclose (gd) + else { + if (gd != NULL) + call gclose (gd) + if (id != NULL) + call gclose (id) + } + if (mgd != NULL) + call gclose (mgd) + if (pfd != NULL) + call close (pfd) + + # Free centering structure. + call apcfree (ap) + + # Close up the file lists. + call imtclose (imlist) + call clpcls (clist) + call clpcls (olist) + call sfree (sp) +end |