diff options
Diffstat (limited to 'pkg/images/tv/imexamine')
31 files changed, 7902 insertions, 0 deletions
diff --git a/pkg/images/tv/imexamine/iecimexam.x b/pkg/images/tv/imexamine/iecimexam.x new file mode 100644 index 00000000..1bcc6d65 --- /dev/null +++ b/pkg/images/tv/imexamine/iecimexam.x @@ -0,0 +1,81 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <error.h> +include "imexam.h" + +# IE_CIMEXAM -- Column plot +# If the input column is INDEF use the last column. + +procedure ie_cimexam (gp, mode, ie, x) + +pointer gp # GIO pointer +int mode # Mode +pointer ie # Structure pointer +real x # Column + +real xavg, junk +int i, x1, x2, y1, y2, nx, ny, npts +pointer sp, title, im, data, ptr, xp, yp + +real asumr() +int clgpseti() +pointer clopset(), ie_gimage(), ie_gdata() +errchk clcpset, clopset + +begin + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + if (IE_PP(ie) != NULL) + call clcpset (IE_PP(ie)) + IE_PP(ie) = clopset ("cimexam") + + if (!IS_INDEF(x)) + IE_X1(ie) = x + + nx = clgpseti (IE_PP(ie), "naverage") + x1 = IE_X1(ie) - (nx - 1) / 2 + 0.5 + x2 = IE_X1(ie) + nx / 2 + 0.5 + xavg = (x1 + x2) / 2. + y1 = INDEFI + y2 = INDEFI + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + return + } + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + npts = nx * ny + + call smark (sp) + call salloc (title, IE_SZTITLE, TY_CHAR) + call salloc (xp, ny, TY_REAL) + + do i = 1, ny + call ie_mwctran (ie, xavg, real(i), junk, Memr[xp+i-1]) + + if (nx > 1) { + ptr = data + call salloc (yp, ny, TY_REAL) + do i = 1, ny { + Memr[yp+i-1] = asumr (Memr[ptr], nx) + ptr = ptr + nx + } + call adivkr (Memr[yp], real (nx), Memr[yp], ny) + } else + yp = data + + call sprintf (Memc[title], IE_SZTITLE, "%s: Columns %d - %d\n%s") + call pargstr (IE_IMNAME(ie)) + call pargi (x1) + call pargi (x2) + call pargstr (IM_TITLE(im)) + + call ie_graph (gp, mode, IE_PP(ie), Memc[title], Memr[xp], + Memr[yp], ny, IE_YLABEL(ie), IE_YFORMAT(ie)) + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/iecolon.x b/pkg/images/tv/imexamine/iecolon.x new file mode 100644 index 00000000..72925500 --- /dev/null +++ b/pkg/images/tv/imexamine/iecolon.x @@ -0,0 +1,1038 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <error.h> +include "imexam.h" + +# List of boundary types, marker types, and colon commands. + +define BTYPES "|constant|nearest|reflect|wrap|project|" +define MTYPES "|point|box|plus|cross|circle|hebar|vebar|hline|vline|diamond|" +define CMDS "|angh|angv|background|banner|boundary|box|buffer|ceiling|\ + |center|constant|dashpat|defkey|eparam|fill|floor|interval|\ + |label|logfile|logx|logy|magzero|majrx|majry|marker|minrx|\ + |minry|naverage|ncolumns|ncontours|ncstat|nhi|nlines|nlstat|\ + |pointmode|radius|round|rplot|select|szmarker|ticklabels|\ + |title|width|x|xlabel|xorder|y|ylabel|yorder|zero|unlearn|\ + |autoredraw|nbins|z1|z2|autoscale|top_closed|allframes|wcs|\ + |xformat|yformat|fitplot|sigma|axes|fittype|beta|iterations|\ + |output|ncoutput|nloutput|" + +define ANGH 1 +define ANGV 2 +define BACKGROUND 3 +define BANNER 4 +define BOUNDARY 5 +define BOX 6 +define BUFFER 7 +define CEILING 8 + +define CENTER 10 +define CONSTANT 11 +define DASHPAT 12 +define DEFKEY 13 +define EPARAM 14 +define FILL 15 +define FLOOR 16 +define INTERVAL 17 + +define LABEL 19 +define LOGFILE 20 +define LOGX 21 +define LOGY 22 +define MAGZERO 23 +define MAJRX 24 +define MAJRY 25 +define MARKER 26 +define MINRX 27 + +define MINRY 29 +define NAVERAGE 30 +define NCOLUMNS 31 +define NCONTOURS 32 +define NCSTAT 33 +define NHI 34 +define NLINES 35 +define NLSTAT 36 + +define POINTMODE 38 +define RADIUS 39 +define ROUND 40 +define RPLOT 41 +define SELECT 42 +define SZMARKER 43 +define TICKLABELS 44 + +define TITLE 46 +define WIDTH 47 +define X 48 +define XLABEL 49 +define XORDER 50 +define Y 51 +define YLABEL 52 +define YORDER 53 +define ZERO 54 +define UNLEARN 55 + +define AUTOREDRAW 57 +define NBINS 58 +define Z1 59 +define Z2 60 +define AUTOSCALE 61 +define TOP_CLOSED 62 +define ALLFRAMES 63 +define WCS 64 + +define XFORMAT 66 +define YFORMAT 67 +define FITPLOT 68 +define SIGMA 69 +define AXES 70 +define FITTYPE 71 +define BETA 72 +define ITERATIONS 73 + +define OUTPUT 75 +define NCOUTPUT 76 +define NLOUTPUT 77 + + +# IE_COLON -- Respond to colon commands. + +procedure ie_colon (ie, cmdstr, gp, redraw) + +pointer ie # IMEXAM data structure +char cmdstr[ARB] # Colon command +pointer gp # GIO pointer +int redraw # Redraw graph? + +char gtype +bool bval +real rval1 +int ival, ncmd +pointer sp, cmd, pp + +bool clgetb(), clgpsetb() +char clgetc() +real clgetr(), clgpsetr() +int nscan(), strdic(), clgeti() +pointer clopset() +errchk clopset, clppsetb, clppsetr, clputb, clputi, clputr + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + # Scan the command string and get the first word. + call sscan (cmdstr) + call gargwrd (Memc[cmd], SZ_LINE) + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, CMDS) + if (ncmd == 0) { + call printf ("Unrecognized or ambiguous command\007") + call sfree (sp) + return + } + + gtype = IE_GTYPE(ie) + pp = IE_PP(ie) + + # Special optimization for the a key. + switch (ncmd) { + case BACKGROUND, CENTER, NAVERAGE, RPLOT, XORDER, WIDTH: + if (IE_LASTKEY(ie) == 'a') { + gtype = 'r' + pp = clopset ("rimexam") + } + if (IE_LASTKEY(ie) == ',') { + gtype = '.' + pp = clopset ("rimexam") + } + } + + # Switch on the command and possibly read further arguments. + switch (ncmd) { + case ANGH: + call gargr (rval1) + if (nscan() == 1) { + call printf ("angh %g\n") + call pargr (clgetr ("simexam.angh")) + } else { + call clputr ("simexam.angh", rval1) + if (gtype == 's') + redraw = YES + } + case ANGV: + call gargr (rval1) + if (nscan() == 1) { + call printf ("angv %g\n") + call pargr (clgetr ("simexam.angv")) + } else { + call clputr ("simexam.angv", rval1) + if (gtype == 's') + redraw = YES + } + case BACKGROUND: + switch (gtype) { + case 'j', 'k', 'r', '.': + call gargb (bval) + if (nscan() == 1) { + call printf ("background %b\n") + call pargb (clgpsetb (pp, "background")) + } else { + call clppsetb (pp, "background", bval) + if (pp == IE_PP(ie)) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case BANNER: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + call gargb (bval) + if (nscan() == 2) { + call clppsetb (pp, "banner", bval) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case BOUNDARY: + call gargwrd (Memc[cmd], SZ_LINE) + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, BTYPES) + if (ncmd == 0) { + call printf ("Boundary types are %s\n") + call pargstr (BTYPES) + } else + call clpstr ("vimexam.boundary", Memc[cmd]) + case BOX: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + call gargb (bval) + if (nscan() == 2) { + call clppsetb (pp, "box", bval) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case BUFFER: + call gargr (rval1) + if (nscan() == 1) { + call printf ("buffer %g\n") + call pargr (clgetr ("rimexam.buffer")) + } else { + call clputr ("rimexam.buffer", rval1) + if (gtype == 'r' || gtype == '.') + redraw = YES + } + case CEILING: + switch (gtype) { + case 's', 'e': + call gargr (rval1) + if (nscan() == 1) { + call printf ("ceiling %g\n") + call pargr (clgpsetr (pp, "ceiling")) + } else { + call clppsetr (pp, "ceiling", rval1) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case CENTER: + switch (gtype) { + case 'j', 'k', 'r', '.': + call gargb (bval) + if (nscan() == 1) { + call printf ("center %b\n") + call pargb (clgpsetb (pp, "center")) + } else { + call clppsetb (pp, "center", bval) + if (pp == IE_PP(ie)) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case CONSTANT: + call gargr (rval1) + if (nscan() == 1) { + call printf ("constant %g\n") + call pargr (clgetr ("vimexam.constant")) + } else + call clputr ("vimexam.constant", rval1) + case DASHPAT: + call gargi (ival) + if (nscan() == 1) { + call printf ("dashpat %g\n") + call pargi (clgeti ("eimexam.dashpat")) + } else { + call clputi ("eimexam.dashpat", ival) + if (gtype == 'e') + redraw = YES + } + case DEFKEY: + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) { + call printf ("defkey %c\n") + call pargc (clgetc ("defkey")) + } else + call clputc ("defkey", Memc[cmd]) + case EPARAM: + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) + Memc[cmd] = gtype + + switch (Memc[cmd]) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 's', 'h', '.': + call gdeactivate (gp, 0) + switch (Memc[cmd]) { + case 'c': + call clcmdw ("eparam cimexam") + case 'j': + call clcmdw ("eparam jimexam") + case 'k': + call clcmdw ("eparam kimexam") + case 'l': + call clcmdw ("eparam limexam") + case 'r', '.': + call clcmdw ("eparam rimexam") + case 's': + call clcmdw ("eparam simexam") + case 'u', 'v': + call clcmdw ("eparam vimexam") + case 'e': + call clcmdw ("eparam eimexam") + case 'h': + call clcmdw ("eparam himexam") + } + if (Memc[cmd] == gtype) + redraw = YES + } + case FILL: + call gargb (bval) + if (nscan() == 1) { + call printf ("fill %b\n") + call pargb (clgetb ("eimexam.fill")) + } else { + call clputb ("eimexam.fill", bval) + if (gtype == 'e') + redraw = YES + } + case FLOOR: + switch (gtype) { + case 's', 'e': + call gargr (rval1) + if (nscan() == 1) { + call printf ("floor %g\n") + call pargr (clgpsetr (pp, "floor")) + } else { + call clppsetr (pp, "floor", rval1) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case INTERVAL: + call gargr (rval1) + if (nscan() == 1) { + call printf ("interval %g\n") + call pargr (clgetr ("eimexam.interval")) + } else { + call clputr ("eimexam.interval", rval1) + if (gtype == 'e') + redraw = YES + } + case LABEL: + call gargb (bval) + if (nscan() == 2) { + call clputb ("eimexam.label", bval) + if (gtype == 'e') + redraw = YES + } + + case LOGFILE: + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) { + call strcpy (IE_LOGFILE(ie), Memc[cmd], SZ_LINE) + if (IE_LOGFD(ie) == NULL) { + call printf ("logfile %s [closed]\n") + call pargstr (Memc[cmd]) + } else { + call printf ("logfile %s [open]\n") + call pargstr (Memc[cmd]) + } + } else { + call clpstr ("logfile", Memc[cmd]) + if (IE_LOGFD(ie) != NULL) { + call close (IE_LOGFD(ie)) + IE_LOGFD(ie) = NULL + } + + call clgstr ("logfile", IE_LOGFILE(ie), SZ_LINE) + if (clgetb ("keeplog")) + iferr (call ie_openlog (ie)) + call erract (EA_WARN) + } + + case LOGX: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'h', '.': + call gargb (bval) + if (nscan() == 2) { + call clppsetb (pp, "logx", bval) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case LOGY: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'h', '.': + call gargb (bval) + if (nscan() == 2) { + call clppsetb (pp, "logy", bval) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case MAGZERO: + call gargr (rval1) + if (nscan() == 1) { + call printf ("magzero %g\n") + call pargr (clgetr ("rimexam.magzero")) + } else { + call clputr ("rimexam.magzero", rval1) + if (gtype == 'r' || gtype == '.') + redraw = YES + } + case AUTOREDRAW: + call gargb (bval) + if (nscan() == 1) { + call printf ("autoredraw %b\n") + call pargb (clgetb ("autoredraw")) + } else + call clputb ("autoredraw", bval) + default: + call ie_colon1 (ie, ncmd, gp, pp, gtype, redraw) + } + + if (pp != IE_PP(ie)) + call clcpset (pp) + if (redraw == YES && !clgetb ("autoredraw")) + redraw = NO + call sfree (sp) +end + + +# IE_COLON1 -- Subprocedure to get around too many strings error in xc. + +procedure ie_colon1 (ie, ncmd, gp, pp, gtype, redraw) + +pointer ie # IMEXAM data structure +int ncmd # Command number +pointer gp # GIO pointer +pointer pp # Pset pointer +char gtype # Graph type +int redraw # Redraw graph? + +int ival +real rval1, rval2 +bool bval +pointer sp, cmd, im + +real clgetr(), clgpsetr() +pointer ie_gimage() +int nscan(), strdic(), clgeti(), clgpseti() +errchk ie_gimage, clppseti + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + switch (ncmd) { + case MAJRX: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + call gargi (ival) + if (nscan() == 1) { + call printf ("majrx %d\n") + call pargi (clgpseti (pp, "majrx")) + } else { + call clppseti (pp, "majrx", ival) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case MAJRY: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + call gargi (ival) + if (nscan() == 1) { + call printf ("majry %d\n") + call pargi (clgpseti (pp, "majry")) + } else { + call clppseti (pp, "majry", ival) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case MARKER: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'h', '.': + call gargwrd (Memc[cmd], SZ_LINE) + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, MTYPES) + if (ncmd == 0) { + call printf ("Marker types are %s\n") + call pargstr (MTYPES) + } else { + call clppset (pp, "marker", Memc[cmd]) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case MINRX: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + call gargi (ival) + if (nscan() == 1) { + call printf ("minrx %d\n") + call pargi (clgpseti (pp, "minrx")) + } else { + call clppseti (pp, "minrx", ival) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case MINRY: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + call gargi (ival) + if (nscan() == 1) { + call printf ("minry %d\n") + call pargi (clgpseti (pp, "minry")) + } else { + call clppseti (pp, "minry", ival) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case NAVERAGE: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'v': + call gargi (ival) + if (nscan() == 1) { + call printf ("naverage %d\n") + call pargi (clgpseti (pp, "naverage")) + } else { + call clppseti (pp, "naverage", ival) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case NCOLUMNS: + switch (gtype) { + case 's', 'e', 'h': + call gargi (ival) + if (nscan() == 1) { + call printf ("ncolumns %d\n") + call pargi (clgpseti (pp, "ncolumns")) + } else { + call clppseti (pp, "ncolumns", ival) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case NCONTOURS: + call gargi (ival) + if (nscan() == 1) { + call printf ("ncontours %g\n") + call pargi (clgeti ("eimexam.ncontours")) + } else { + call clputi ("eimexam.ncontours", ival) + if (gtype == 'e') + redraw = YES + } + case NCSTAT: + call gargi (ival) + if (nscan() == 1) { + call printf ("ncstat %g\n") + call pargi (clgeti ("ncstat")) + } else + call clputi ("ncstat", ival) + case NHI: + call gargi (ival) + if (nscan() == 1) { + call printf ("nhi %g\n") + call pargi (clgeti ("eimexam.nhi")) + } else { + call clputi ("eimexam.nhi", ival) + if (gtype == 'e') + redraw = YES + } + case NLINES: + switch (gtype) { + case 's', 'e', 'h': + call gargi (ival) + if (nscan() == 1) { + call printf ("nlines %d\n") + call pargi (clgpseti (pp, "nlines")) + } else { + call clppseti (pp, "nlines", ival) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case NLSTAT: + call gargi (ival) + if (nscan() == 1) { + call printf ("nlstat %g\n") + call pargi (clgeti ("nlstat")) + } else + call clputi ("nlstat", ival) + case POINTMODE: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'h', '.': + call gargb (bval) + if (nscan() == 2) { + call clppsetb (pp, "pointmode", bval) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case RADIUS: + call gargr (rval1) + if (nscan() == 1) { + call printf ("radius %g\n") + call pargr (clgetr ("rimexam.radius")) + } else { + call clputr ("rimexam.radius", rval1) + if (gtype == 'r' || gtype == '.') + redraw = YES + } + case ROUND: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + call gargb (bval) + if (nscan() == 2) { + call clppsetb (pp, "round", bval) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case RPLOT: + switch (gtype) { + case 'j', 'k', 'r', '.': + call gargr (rval1) + if (nscan() == 1) { + call printf ("rplot %g\n") + call pargr (clgpsetr (pp, "rplot")) + } else { + call clppsetr (pp, "rplot", rval1) + if (pp == IE_PP(ie)) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case SELECT: + call gargi (ival) + if (nscan () > 1) { + if (IE_LIST(ie) != NULL) + IE_INDEX(ie) = ival + else + IE_NEWFRAME(ie) = ival + IE_MAPFRAME(ie) = 0 + iferr (im = ie_gimage (ie, YES)) + call erract (EA_WARN) + } + case SZMARKER: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'h', '.': + call gargi (ival) + if (nscan() == 1) { + call printf ("szmarker %d\n") + call pargi (clgpseti (pp, "szmarker")) + } else { + call clppseti (pp, "szmarker", ival) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case TICKLABELS: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + call gargb (bval) + if (nscan() == 2) { + call clppsetb (pp, "ticklabels", bval) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case TITLE: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 's', 'v', 'e', 'h', '.': + Memc[cmd] = EOS + call gargstr (Memc[cmd], SZ_LINE) + call clppset (pp, "title", Memc[cmd]) + redraw = YES + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case WIDTH: + switch (gtype) { + case 'j', 'k', 'r', '.': + call gargr (rval1) + if (nscan() == 1) { + call printf ("width %g\n") + call pargr (clgpsetr (pp, "width")) + } else { + call clppsetr (pp, "width", rval1) + if (pp == IE_PP(ie)) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case X: + switch (gtype) { + case 'c', 'j', 'k', 'l', 'r', 'v', 'h', '.': + call gargr (rval1) + call gargr (rval2) + if (nscan() < 3) { + call clppsetr (pp, "x1", INDEF) + call clppsetr (pp, "x2", INDEF) + } else { + call clppsetr (pp, "x1", rval1) + call clppsetr (pp, "x2", rval2) + } + redraw = YES + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case XLABEL: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + Memc[cmd] = EOS + call gargstr (Memc[cmd], SZ_LINE) + call clppset (pp, "xlabel", Memc[cmd]) + redraw = YES + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case XORDER: + switch (gtype) { + case 'j', 'k', 'r', '.': + call gargi (ival) + if (nscan() == 1) { + call printf ("xorder %d\n") + call pargi (clgpseti (pp, "xorder")) + } else { + call clppseti (pp, "xorder", ival) + if (pp == IE_PP(ie)) + redraw = YES + } + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case Y: + switch (gtype) { + case 'c', 'j', 'k', 'l', 'r', 'v', 'h', '.': + call gargr (rval1) + call gargr (rval2) + if (nscan() < 3) { + call clppsetr (pp, "y1", INDEF) + call clppsetr (pp, "y2", INDEF) + } else { + call clppsetr (pp, "y1", rval1) + call clppsetr (pp, "y2", rval2) + } + redraw = YES + default: + call printf ("Parameter does not apply to current graph\007\n") + } + default: + call ie_colon2 (ie, ncmd, gp, pp, gtype, redraw) + } + + call sfree (sp) +end + + +# IE_COLON2 -- Subprocedure to get around too many strings error in xc. + +procedure ie_colon2 (ie, ncmd, gp, pp, gtype, redraw) + +pointer ie # IMEXAM data structure +int ncmd # Command number +pointer gp # GIO pointer +pointer pp # Pset pointer +char gtype # Graph type +int redraw # Redraw graph? + +int ival +real rval1 +bool bval +pointer sp, cmd + +real clgetr() +bool clgetb() +int nscan(), clgeti(), btoi(), strdic() + +begin + call smark (sp) + call salloc (cmd, SZ_LINE, TY_CHAR) + + switch (ncmd) { + case YLABEL: + switch (gtype) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 'h', '.': + Memc[cmd] = EOS + call gargstr (Memc[cmd], SZ_LINE) + call clppset (pp, "ylabel", Memc[cmd]) + redraw = YES + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case YORDER: + call gargi (ival) + if (nscan() == 1) { + call printf ("yorder %d\n") + call pargi (clgeti ("rimexam.yorder")) + } else { + call clputi ("rimexam.yorder", ival) + if (gtype == 'r' || gtype == '.') + redraw = YES + } + case ZERO: + call gargr (rval1) + if (nscan() == 1) { + call printf ("zero %g\n") + call pargr (clgetr ("eimexam.zero")) + } else { + call clputr ("eimexam.zero", rval1) + if (gtype == 'e') + redraw = YES + } + case UNLEARN: + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) + Memc[cmd] = gtype + + switch (Memc[cmd]) { + case 'c', 'u', 'j', 'k', 'l', 'r', 'v', 'e', 's', 'h', '.': + switch (Memc[cmd]) { + case 'c': + call clcmdw ("unlearn cimexam") + case 'j': + call clcmdw ("unlearn jimexam") + case 'k': + call clcmdw ("unlearn jimexam") + case 'l': + call clcmdw ("unlearn limexam") + case 'r', '.': + call clcmdw ("unlearn rimexam") + case 's': + call clcmdw ("unlearn simexam") + case 'u', 'v': + call clcmdw ("unlearn vimexam") + case 'e': + call clcmdw ("unlearn eimexam") + case 'h': + call clcmdw ("unlearn himexam") + } + if (Memc[cmd] == gtype) + redraw = YES + default: + call printf ("Parameter does not apply to current graph\007\n") + } + case NBINS: + call gargi (ival) + if (nscan() == 1) { + call printf ("nbins %d\n") + call pargi (clgeti ("himexam.nbins")) + } else { + call clputi ("himexam.nbins", ival) + if (gtype == 'h') + redraw = YES + } + case Z1: + call gargr (rval1) + if (nscan() == 1) { + call printf ("z1 %g\n") + call pargr (clgetr ("himexam.z1")) + } else { + call clputr ("himexam.z1", rval1) + if (gtype == 'h') + redraw = YES + } + case Z2: + call gargr (rval1) + if (nscan() == 1) { + call printf ("z2 %g\n") + call pargr (clgetr ("himexam.z2")) + } else { + call clputr ("himexam.z2", rval1) + if (gtype == 'h') + redraw = YES + } + case AUTOSCALE: + call gargb (bval) + if (nscan() == 1) { + call printf ("autoscale %b\n") + call pargb (clgetb ("himexam.autoscale")) + } else { + call clputb ("himexam.autoscale", bval) + if (gtype == 'h') + redraw = YES + } + case TOP_CLOSED: + call gargb (bval) + if (nscan() == 1) { + call printf ("top_closed %b\n") + call pargb (clgetb ("himexam.top_closed")) + } else { + call clputb ("himexam.top_closed", bval) + if (gtype == 'h') + redraw = YES + } + case ALLFRAMES: + call gargb (bval) + if (nscan() == 1) { + call printf ("allframes %b\n") + call pargb (clgetb ("allframes")) + } else { + call clputb ("allframes", bval) + IE_ALLFRAMES(ie) = btoi (bval) + } + case WCS: + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) { + call printf ("wcs %s\n") + call pargstr (IE_WCSNAME(ie)) + } else { + call strcpy (Memc[cmd], IE_WCSNAME(ie), SZ_FNAME) + call ie_mwinit (ie) + redraw = YES + } + case XFORMAT: + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) + call clpstr ("xformat", "") + else + call clpstr ("xformat", Memc[cmd]) + case YFORMAT: + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) + call clpstr ("yformat", "") + else + call clpstr ("yformat", Memc[cmd]) + case FITPLOT: + call gargb (bval) + if (nscan() == 1) { + call printf ("fitplot %b\n") + call pargb (clgetb ("rimexam.fitplot")) + } else { + call clputb ("rimexam.fitplot", bval) + if (gtype == 'r') + redraw = YES + } + case SIGMA: + call gargr (rval1) + if (nscan() == 1) { + call printf ("sigma %g\n") + call pargr (clgetr ("jimexam.sigma")) + } else { + call clputr ("jimexam.sigma", rval1) + if (gtype == 'j' || gtype == 'k') + redraw = YES + } + case AXES: + call gargb (bval) + if (nscan() == 2) { + call clputb ("simexam.axes", bval) + if (gtype == 's') + redraw = YES + } + case FITTYPE: + call gargwrd (Memc[cmd], SZ_LINE) + if (nscan() == 1) { + call clgstr ("rimexam.fittype", Memc[cmd], SZ_LINE) + call printf ("fittype %s\n") + call pargstr (Memc[cmd]) + } else { + ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, + "|gaussian|moffat|") + if (ncmd == 0) { + call printf ("Profile fit types are %s\n") + call pargstr ("|gaussian|moffat|") + } else { + call clpstr ("rimexam.fittype", Memc[cmd]) + if (gtype == 'r' || gtype == '.') + redraw = YES + } + } + case BETA: + call gargr (rval1) + if (nscan() == 1) { + call printf ("beta %g\n") + call pargr (clgetr ("rimexam.beta")) + } else { + call clputr ("rimexam.beta", rval1) + if (gtype == 'r' || gtype == '.') + redraw = YES + } + case ITERATIONS: + call gargi (ival) + if (nscan() == 1) { + call printf ("iterations %d\n") + call pargi (clgeti ("rimexam.iterations")) + } else { + call clputi ("rimexam.iterations", ival) + if (gtype == 'r') + redraw = YES + } + + case OUTPUT: + call gargwrd (Memc[cmd], SZ_FNAME) + if (nscan() == 1) { + call clgstr ("output", Memc[cmd], SZ_FNAME) + call printf ("output `%s'\n") + call pargstr (Memc[cmd]) + } else + call clpstr ("output", Memc[cmd]) + case NCOUTPUT: + call gargi (ival) + if (nscan() == 1) { + call printf ("ncoutput %g\n") + call pargi (clgeti ("ncoutput")) + } else + call clputi ("ncoutput", ival) + case NLOUTPUT: + call gargi (ival) + if (nscan() == 1) { + call printf ("nloutput %g\n") + call pargi (clgeti ("nloutput")) + } else + call clputi ("nloutput", ival) + + default: + call printf ("Ambiguous or unrecognized command\007\n") + } + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/iedisplay.x b/pkg/images/tv/imexamine/iedisplay.x new file mode 100644 index 00000000..4015bca7 --- /dev/null +++ b/pkg/images/tv/imexamine/iedisplay.x @@ -0,0 +1,55 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> + +# IE_DISPLAY -- Display an image. For the sake of convenience in this +# prototype program we do this by calling a task via the cl. This is an +# interface violation which we try to mitigate by using a CL parameter to +# hide the knowledge of how to format the command (as well as make it easy +# for the user to control how images are displayed). + +procedure ie_display (ie, image, frame) + +pointer ie #I imexamine descriptor +char image[ARB] #I image to be displayed +int frame #I frame in which to display image + +int nchars +pointer sp, d_cmd, d_args, d_template, im +int gstrcpy(), strmac(), ie_getnframes() +pointer immap() + +begin + call smark (sp) + call salloc (d_cmd, SZ_LINE, TY_CHAR) + call salloc (d_args, SZ_LINE, TY_CHAR) + call salloc (d_template, SZ_LINE, TY_CHAR) + + # Verify that the named image or image section exists. + iferr (im = immap (image, READ_ONLY, 0)) { + call erract (EA_WARN) + call sfree (sp) + return + } else + call imunmap (im) + + # Get the display command template. + call clgstr ("display", Memc[d_template], SZ_LINE) + + # Construct the macro argument list, a sequence of EOS delimited + # strings terminated by a double EOS. + + call aclrc (Memc[d_args], SZ_LINE) + nchars = gstrcpy (image, Memc[d_args], SZ_LINE) + 1 + call sprintf (Memc[d_args+nchars], SZ_LINE-nchars, "%d") + call pargi (frame) + + # Expand the command template to form the CL command. + nchars = strmac (Memc[d_template], Memc[d_args], Memc[d_cmd], SZ_LINE) + + # Send the command off to the CL and wait for completion. + call clcmdw (Memc[d_cmd]) + nchars = ie_getnframes (ie) + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/ieeimexam.x b/pkg/images/tv/imexamine/ieeimexam.x new file mode 100644 index 00000000..059721ba --- /dev/null +++ b/pkg/images/tv/imexamine/ieeimexam.x @@ -0,0 +1,243 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <gset.h> +include <config.h> +include <mach.h> +include <imhdr.h> +include <xwhen.h> +include <fset.h> +include "imexam.h" + + +# IE_EIMEXAM -- Contour map +# This is an interface to the NCAR CONREC routine. + +procedure ie_eimexam (gp, mode, ie, x, y) + +pointer gp # GIO pointer +int mode # Mode +pointer ie # IE pointer +real x, y # Center + +bool banner +int nset, ncontours, dashpat, nhi +int x1, x2, y1, y2, nx, ny, npts, wkid +real vx1, vx2, vy1, vy2, xs, xe, ys, ye +real interval, floor, ceiling, zero, finc, zmin, zmax +pointer sp, title, hostid, user, xlabel, ylabel, im, data, data1 + +pointer pp, clopset(), ie_gdata(), ie_gimage() +bool clgpsetb(), fp_equalr() +int clgpseti(), btoi() +real clgpsetr() + +int isizel, isizem, isizep, nrep, ncrt, ilab, nulbll, ioffd +int ioffm, isolid, nla, nlm +real xlt, ybt, side, ext, hold[5] +common /conre4/ isizel, isizem , isizep, nrep, ncrt, ilab, nulbll, + ioffd, ext, ioffm, isolid, nla, nlm, xlt, ybt, side +int first +common /conflg/ first +common /noaolb/ hold + +begin + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + pp = IE_PP(ie) + if (pp != NULL) + call clcpset (pp) + pp = clopset ("eimexam") + IE_PP(ie) = pp + + if (!IS_INDEF(x)) + IE_X1(ie) = x + if (!IS_INDEF(y)) + IE_Y1(ie) = y + + nx = clgpseti (pp, "ncolumns") + ny = clgpseti (pp, "nlines") + x1 = IE_X1(ie) - (nx - 1) / 2 + 0.5 + x2 = IE_X1(ie) + nx / 2 + 0.5 + y1 = IE_Y1(ie) - (ny - 1) / 2 + 0.5 + y2 = IE_Y1(ie) + ny / 2 + 0.5 + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + return + } + + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + npts = nx * ny + xs = x1 + xe = x2 + ys = y1 + ye = y2 + + call smark (sp) + banner = false + if (mode == NEW_FILE) { + call gclear (gp) + + # Set the WCS + call gswind (gp, xs, xe, ys, ye) + + if (!clgpsetb (pp, "fill")) + call gsetr (gp, G_ASPECT, real (ny-1) / real (nx-1)) + + call gseti (gp, G_ROUND, btoi (clgpsetb (pp, "round"))) + + if (clgpsetb (pp, "box")) { + # Get number of major and minor tick marks. + call gseti (gp, G_XNMAJOR, clgpseti (pp, "majrx")) + call gseti (gp, G_XNMINOR, clgpseti (pp, "minrx")) + call gseti (gp, G_YNMAJOR, clgpseti (pp, "majry")) + call gseti (gp, G_YNMINOR, clgpseti (pp, "minry")) + + # Label tick marks on axes? + call gseti (gp, G_LABELTICKS, + btoi (clgpsetb (pp, "ticklabels"))) + + # Labels + call salloc (title, IE_SZTITLE, TY_CHAR) + call salloc (hostid, SZ_LINE, TY_CHAR) + call salloc (user, SZ_LINE, TY_CHAR) + call salloc (xlabel, SZ_LINE, TY_CHAR) + call salloc (ylabel, SZ_LINE, TY_CHAR) + + banner = clgpsetb (pp, "banner") + if (banner) { + call sysid (Memc[hostid], SZ_LINE) + # We must postpone the parameter line until after conrec. + call sprintf (Memc[title], IE_SZTITLE, "%s\n\n%s") + call pargstr (Memc[hostid]) + call pargstr (IM_TITLE(im)) + } else + Memc[title] = EOS + + call clgpset (pp, "title", Memc[user], SZ_LINE) + if (Memc[user] != EOS) { + call strcat ("\n", Memc[title], IE_SZTITLE) + call strcat (Memc[user], Memc[title], IE_SZTITLE) + } + call clgpset (pp, "xlabel", Memc[xlabel], SZ_LINE) + call clgpset (pp, "ylabel", Memc[ylabel], SZ_LINE) + + call glabax (gp, Memc[title], Memc[xlabel], Memc[ylabel]) + } + } + + # First of all, intialize conrec's block data before altering any + # parameters in common. + first = 1 + call conbd + + # Set contour parameters + zero = clgpsetr (pp, "zero") + floor = clgpsetr (pp, "floor") + ceiling = clgpsetr (pp, "ceiling") + nhi = clgpseti (pp, "nhi") + dashpat = clgpseti (pp, "dashpat") + + # Resolve INDEF limits. + if (IS_INDEF (floor) || IS_INDEF (ceiling)) { + call alimr (Memr[data], npts, zmin, zmax) + if (IS_INDEF (floor)) + floor = zmin + if (IS_INDEF (ceiling)) + ceiling = zmax + } + + # Apply the zero point shift. + if (abs (zero) > EPSILON) { + call salloc (data1, npts, TY_REAL) + call asubkr (Memr[data], zero, Memr[data1], npts) + floor = floor - zero + ceiling = ceiling - zero + } else + data1 = data + + # Avoid conrec's automatic scaling. + if (floor == 0.) + floor = EPSILON + if (ceiling == 0.) + ceiling = EPSILON + + # The user can suppress the contour labelling by setting the common + # parameter "ilab" to zero. + if (btoi (clgpsetb (pp, "label")) == NO) + ilab = 0 + else + ilab = 1 + + # User can specify either the number of contours or the contour + # interval, or let conrec pick a nice number. Get params and + # encode the FINC param expected by conrec. + + ncontours = clgpseti (pp, "ncontours") + if (ncontours <= 0) { + interval = clgpsetr (pp, "interval") + if (interval <= 0) + finc = 0 + else + finc = interval + } else + finc = - abs (ncontours) + + # Open device and make contour plot. + call gopks (STDERR) + wkid = 1 + call gopwk (wkid, 6, gp) + call gacwk (wkid) + + # Make the contour plot. + nset = 1 # No conrec viewport + ioffm = 1 # No conrec box + call gswind (gp, 1., real (nx), 1., real (ny)) + call ggview (gp, vx1, vx2, vy1, vy2) + call set (vx1, vx2, vy1, vy2, 1.0, real (nx), 1.0, real (ny), 1) + call conrec (Memr[data1], nx, nx, ny, floor, + ceiling, finc, nset, nhi, -dashpat) + + call gdawk (wkid) + call gclks () + + call gswind (gp, xs, xe, ys, ye) + if (banner) { + if (fp_equalr (hold(5), 1.0)) { + call sprintf (Memc[title], IE_SZTITLE, + "%s\n%s: Contoured from %g to %g, interval = %g\n%s") + call pargstr (Memc[hostid]) + call pargstr (IE_IMNAME(ie)) + call pargr (hold(1)) + call pargr (hold(2)) + call pargr (hold(3)) + call pargstr (IM_TITLE(im)) + } else { + call sprintf (Memc[title], IE_SZTITLE, + "%s\n%s:contoured from %g to %g, interval = %g, labels scaled by %g\n%s") + call pargstr (Memc[xlabel]) + call pargstr (IE_IMNAME(ie)) + call pargr (hold(1)) + call pargr (hold(2)) + call pargr (hold(3)) + call pargr (hold(5)) + call pargstr (IM_TITLE(im)) + } + + if (Memc[user] != EOS) { + call strcat ("\n", Memc[user], IE_SZTITLE) + call strcat (Memc[user], Memc[title], IE_SZTITLE) + } + + call gseti (gp, G_DRAWAXES, NO) + call glabax (gp, Memc[title], "", "") + + } else + call gtext (gp, xs, ys, "", "") + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/iegcur.x b/pkg/images/tv/imexamine/iegcur.x new file mode 100644 index 00000000..2b76cee5 --- /dev/null +++ b/pkg/images/tv/imexamine/iegcur.x @@ -0,0 +1,242 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <ctype.h> +include <mach.h> +include "imexam.h" + +# IE_GCUR -- Get IMEXAM cursor value. +# This is an interface between the standard cursor input and IMEXAM. +# It reads the appropriate cursor, determines the image index or frame +# type, makes the appropriate default coordinate conversions when using +# graphics cursor input, and gets any further cursor reads needed. +# Missing coordinates default to the last coordinates. + +int procedure ie_gcur (ie, curtype, x, y, key, strval, maxch) + +pointer ie #I IMEXAM structure +int curtype #I cursor type (0=image, 1=graphics, 2=text) +real x, y #O cursor position +int key #O keystroke value of cursor event +char strval[ARB] #O string value, if any +int maxch #I max chars out + +char ch +real x1, y1, x2, y2, dx, dy, r, cosa, sina +int temp, k[2], nitems, wcs, ip, i + +bool streq() +char clgetc() +int clgcur(), imd_gcur(), ctor(), cctoc() +errchk clgcur, imd_gcur + +begin + # Save last cursor value. + x1 = x; y1 = y + strval[1] = EOS + k[1] = clgetc ("defkey") + + # Get one or more cursor values from the desired cursor parameter. + # Check for missing coordinates and substitute the last value. + + do i = 1, 2 { + switch (curtype) { + case 'i': + nitems = imd_gcur ("imagecur", x, y, wcs, k[i], strval, maxch) + if (IS_INDEF(x)) + x = x1 + if (IS_INDEF(y)) + y = y1 + IE_NEWFRAME(ie) = wcs + if (IE_DFRAME(ie) <= 0) + IE_DFRAME(ie) = IE_NEWFRAME(ie) + + case 'g': + nitems = clgcur ("graphcur", x, y, wcs, k[i], strval, maxch) + + # Make any needed default coordinate conversions from the + # graphic coordinates. + + switch (IE_GTYPE(ie)) { + case 'c', 'k': # Column plot + y = x + x = IE_X1(ie) + + if (IS_INDEF(y)) + y = y1 + else if (IE_MW(ie) != NULL) { + if (streq (IE_WCSNAME(ie), "logical")) + ; + else if (streq (IE_WCSNAME(ie), "physical")) + call ie_imwctran (ie, x, y, dx, y) + else { + r = y + y = IM_LEN(IE_IM(ie),2) + call ie_mwctran (ie, x, 1., dx, y1) + call ie_mwctran (ie, x, y, dx, y2) + dy = y + while (dy > .001) { + dy = dy / 2 + if (r > y1) { + if (r < y2) + y = y - dy + else + y = y + dy + } else { + if (r < y2) + y = y + dy + else + y = y - dy + } + call ie_mwctran (ie, x, y, dx, y2) + } + } + } + case 'e': # Contour plot + if (IS_INDEF(x)) + x = x1 + if (IS_INDEF(y)) + y = y1 + case 'j', 'l': # Line plot + y = IE_Y1(ie) + + if (IS_INDEF(x)) + x = x1 + else if (IE_MW(ie) != NULL) { + if (streq (IE_WCSNAME(ie), "logical")) + ; + else if (streq (IE_WCSNAME(ie), "physical")) + call ie_imwctran (ie, x, y, x, dy) + else { + r = x + x = IM_LEN(IE_IM(ie),1) + call ie_mwctran (ie, 1., y, x1, dy) + call ie_mwctran (ie, x, y, x2, dy) + dx = x + while (dx > .001) { + dx = dx / 2 + if (r > x1) { + if (r < x2) + x = x - dx + else + x = x + dx + } else { + if (r < x2) + x = x + dx + else + x = x - dx + } + call ie_mwctran (ie, x, y, x2, dy) + } + } + } + case 'r','.': # Radial profile plot + x = IE_X1(ie) + y = IE_Y1(ie) + case 'h', 's': # Surface plot + x = IE_X1(ie) + y = IE_Y1(ie) + case 'u': # Vector plot + if (IS_INDEF(x)) + x = x1 + y = x * sina + (IE_Y1(ie) + IE_Y2(ie)) / 2 + x = x * cosa + (IE_X1(ie) + IE_X2(ie)) / 2 + case 'v': # Vector plot + if (IS_INDEF(x)) + x = x1 + y = x * sina + IE_Y1(ie) + x = x * cosa + IE_X1(ie) + } + } + + key = k[1] + switch (key) { + case 'v', 'u': + if (i == 1) { + x1 = x + y1 = y + call printf ("again:") + } else { + x2 = x + y2 = y + r = sqrt (real ((y2-y1)**2 + (x2-x1)**2)) + if (r > 0.) { + cosa = (x2 - x1) / r + sina = (y2 - y1) / r + } else { + cosa = 0. + sina = 0. + } + call printf ("\n") + switch (key) { + case 'v': + x = x1 + y = y1 + case 'u': + x = 2 * x1 - x2 + y = 2 * y1 - y2 + } + IE_X2(ie) = x2 + IE_Y2(ie) = y2 + break + } + case 'b': + if (i == 1) { + IE_IX1(ie) = x + 0.5 + IE_IY1(ie) = y + 0.5 + call printf ("again:") + } else { + IE_IX2(ie) = x + 0.5 + IE_IY2(ie) = y + 0.5 + call printf ("\n") + temp = IE_IX1(ie) + IE_IX1(ie) = min (IE_IX1(ie), IE_IX2(ie)) + IE_IX2(ie) = max (temp, IE_IX2(ie)) + temp = IE_IY1(ie) + IE_IY1(ie) = min (IE_IY1(ie), IE_IY2(ie)) + IE_IY2(ie) = max (temp, IE_IY2(ie)) + break + } + default: + break + } + } + + # Map numeric colon sequences (: x [y] key strval) to make them appear + # as ordinary "x y key" type cursor reads. This makes it possible for + # the user to access any command using typed in rather than positional + # cursor coordinates. Special treatment is also given to the syntax + # ":lN" and ":cN", provided for compatibility with IMPLOT for simple + # line and column plots. + + if (key == ':') { + for (ip=1; IS_WHITE(strval[ip]); ip=ip+1) + ; + if (IS_DIGIT(strval[ip])) { + if (ctor (strval, ip, x) <= 0) + ; + if (ctor (strval, ip, y) <= 0) + y = x + for (; IS_WHITE(strval[ip]); ip=ip+1) + ; + if (cctoc (strval, ip, ch) > 0) + key = ch + call strcpy (strval[ip], strval, maxch) + + } else if (strval[ip] == 'l' && IS_DIGIT(strval[ip+1])) { + ip = ip + 1 + if (ctor (strval, ip, x) > 0) { + y = x + key = 'l' + } + } else if (strval[ip] == 'c' && IS_DIGIT(strval[ip+1])) { + ip = ip + 1 + if (ctor (strval, ip, x) > 0) { + y = x + key = 'c' + } + } + } + + return (nitems) +end diff --git a/pkg/images/tv/imexamine/iegdata.x b/pkg/images/tv/imexamine/iegdata.x new file mode 100644 index 00000000..6e1f7e91 --- /dev/null +++ b/pkg/images/tv/imexamine/iegdata.x @@ -0,0 +1,45 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +# IE_GDATA -- Get image data with boundary checking. + +pointer procedure ie_gdata (im, x1, x2, y1, y2) + +pointer im # IMIO pointer +int x1, x2, y1, y2 # Subraster limits (input and output) + +int i, nc, nl +pointer imgs2r() +errchk imgs2r + +begin + nc = IM_LEN(im,1) + nl = IM_LEN(im,2) + + if (IS_INDEFI (x1)) + x1 = 1 + if (IS_INDEFI (x2)) + x2 = nc + if (IS_INDEFI (y1)) + y1 = 1 + if (IS_INDEFI (y2)) + y2 = nl + + i = max (x1, x2) + x1 = min (x1, x2) + x2 = i + i = max (y1, y2) + y1 = min (y1, y2) + y2 = i + + if (x2 < 1 || x1 > nc || y2 < 1 || y1 > nl) + call error (1, "Pixels out of bounds") + + x1 = max (1, x1) + x2 = min (nc, x2) + y1 = max (1, y1) + y2 = min (nl, y2) + + return (imgs2r (im, x1, x2, y1, y2)) +end diff --git a/pkg/images/tv/imexamine/iegimage.x b/pkg/images/tv/imexamine/iegimage.x new file mode 100644 index 00000000..b0fda919 --- /dev/null +++ b/pkg/images/tv/imexamine/iegimage.x @@ -0,0 +1,261 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include "imexam.h" + +# IE_GIMAGE -- Get input image name and return IMIO pointer. +# If examining a list of images access the indexed image, displaying it if +# not already displayed. Otherwise the image loaded into the current display +# frame is displayed, if it can be accessed, or the image frame buffer itself +# is examined. If there is neither a list of images nor display access the +# user is queried for the name of the image to be examined. +# This procedure uses a prototype display interface (IMD/IW). + +pointer procedure ie_gimage (ie, select) + +pointer ie #I IMEXAM pointer +int select #I select frame? + +char errstr[SZ_FNAME] +int frame, i, j, k +pointer sp, image, dimage, imname, im + +int imtrgetim(), fnldir(), errget() +bool strne(), streq() +pointer imd_mapframe(), immap() +errchk imd_mapframe, immap, ie_display, ie_mwinit + + +begin + call smark (sp) + call salloc (image, SZ_FNAME, TY_CHAR) + call salloc (imname, SZ_FNAME, TY_CHAR) + call salloc (dimage, SZ_FNAME, TY_CHAR) + + # Get image name, and display image if using display. If we are + # examining a list of images, the list and the current index into + # the list determine the image to be examined. If there is no list + # we examine the currently displayed images, if any, else the + # contents of the image display frame buffers are examined as images. + + if (IE_LIST(ie) != NULL) { + # Get image name. + IE_INDEX(ie) = max(1, min(IE_LISTLEN(ie), IE_INDEX(ie))) + if (imtrgetim (IE_LIST(ie), IE_INDEX(ie), Memc[image], + SZ_FNAME) == EOF) + call error (1, "Reference outside of image list") + + # Display image. + if (IE_USEDISPLAY(ie) == YES) { + # Is named image currently loaded into the image display? + frame = 0 + if (streq (Memc[image], IE_IMAGE(ie))) + frame = IE_MAPFRAME(ie) + else { + if (IE_DS(ie) == NULL) + IE_DS(ie) = imd_mapframe (max (1, IE_NEWFRAME(ie)/100), + READ_WRITE, NO) + + do i = 1, IE_NFRAMES(ie) { + if (i == IE_MAPFRAME(ie)/100) + next + do j = 1, 99 { + k = i * 100 + j + iferr (call ie_imname (IE_DS(ie), k, + Memc[dimage], SZ_FNAME)) + break + if (streq (Memc[image], Memc[dimage])) { + frame = k + break + } + } + if (frame != 0) + break + } + } + + # Load image into display frame if not already loaded. + # If the allframes option is specified cycle through the + # available display frames, otherwise resuse the same frame. + + if (frame == 0) { + if (IE_DS(ie) != NULL) { + if (IE_IM(ie) == IE_DS(ie)) + IE_IM(ie) = NULL + call imunmap (IE_DS(ie)) + } + + frame = 100 * max (1, IE_DFRAME(ie) / 100) + 1 + call ie_display (ie, Memc[image], frame/100) + + IE_MAPFRAME(ie) = 0 + if (IE_ALLFRAMES(ie) == YES) { + IE_DFRAME(ie) = frame + 100 + if (IE_DFRAME(ie)/100 > IE_NFRAMES(ie)) + IE_DFRAME(ie) = 101 + } + } + + # Map and display-select the frame. + if (frame != IE_MAPFRAME(ie) || frame != IE_NEWFRAME(ie)) { + if (IE_DS(ie) != NULL) { + if (IE_IM(ie) == IE_DS(ie)) + IE_IM(ie) = NULL + call imunmap (IE_DS(ie)) + } + IE_DS(ie) = imd_mapframe (frame/100, READ_WRITE, select) + IE_MAPFRAME(ie) = frame + IE_NEWFRAME(ie) = frame + } + } + + } else if (IE_USEDISPLAY(ie) == YES) { + # Map the new display frame. + if (IE_NEWFRAME(ie) != IE_MAPFRAME(ie)) { + if (IE_NEWFRAME(ie)/100 != IE_MAPFRAME(ie)/100) { + if (IE_DS(ie) != NULL) { + if (IE_IM(ie) == IE_DS(ie)) + IE_IM(ie) = NULL + call imunmap (IE_DS(ie)) + } + IE_DS(ie) = imd_mapframe (IE_NEWFRAME(ie)/100, READ_WRITE, + select) + } + IE_MAPFRAME(ie) = IE_NEWFRAME(ie) + } + + # Get the image name. + call ie_imname (IE_DS(ie), IE_MAPFRAME(ie), Memc[image], SZ_FNAME) + + } else + call clgstr ("image", Memc[image], SZ_FNAME) + + # Check if the image has not been mapped and if so map it. + # Possibly log any change of image. Always map the physical image, + # not a section, since we do everything in image coordinates. + + if (IE_IM(ie) == NULL || strne (Memc[image], IE_IMAGE(ie))) { + + # Strip the path. + call imgcluster (Memc[image], Memc[imname], SZ_FNAME) + i = fnldir (Memc[imname], Memc[imname], SZ_FNAME) + call strcpy (Memc[image+i], IE_IMNAME(ie), IE_SZFNAME) + + # Map the image. + iferr (im = immap (Memc[image], READ_ONLY, 0)) { + # Warn user once. + i = errget (Memc[imname], SZ_FNAME) + if (strne (Memc[imname], errstr)) { + call erract (EA_WARN) + call strcpy (Memc[imname], errstr, SZ_FNAME) + } + + # Access the display frame buffer as the data image. + if (IE_USEDISPLAY(ie) == YES && IE_LIST(ie) == NULL) { + if (IE_IM(ie) != NULL && IE_IM(ie) != IE_DS(ie)) + iferr (call imunmap (IE_IM(ie))) + ; + IE_IM(ie) = IE_DS(ie) + call sprintf (IE_IMAGE(ie), IE_SZFNAME, "Frame.%d(%s)") + call pargi (IE_MAPFRAME(ie)) + call pargstr (IE_IMNAME(ie)) + call strcpy ("Contents of raw image frame buffer\n", + IM_TITLE(IE_IM(ie)), SZ_IMTITLE) + } else + call erract (EA_WARN) + + } else { + # Adjust image sections. + call ie_gimage1 (im, Memc[image], Memc[imname], SZ_FNAME) + if (strne (Memc[image], Memc[imname])) { + call imunmap (im) + im = immap (Memc[imname], READ_ONLY, 0) + } + + # Make the new image the current one. + errstr[1] = EOS + call strcpy (Memc[image], IE_IMAGE(ie), IE_SZFNAME) + if (IE_IM(ie) != NULL && IE_IM(ie) != IE_DS(ie)) + iferr (call imunmap (IE_IM(ie))) + ; + if (IE_MW(ie) != NULL) + call mw_close (IE_MW(ie)) + IE_IM(ie) = im + if (IE_LOGFD(ie) != NULL) { + call fprintf (IE_LOGFD(ie), "# [%d] %s - %s\n") + call pargi (IE_INDEX(ie)) + call pargstr (IE_IMNAME(ie)) + call pargstr (IM_TITLE(IE_IM(ie))) + } + } + } + + call ie_mwinit (ie) + + call sfree (sp) + return (IE_IM(ie)) +end + + +# IE_GIMAGE1 -- Convert input image section name to a 2D physical image section. + +procedure ie_gimage1 (im, input, output, maxchar) + +pointer im #I IMIO pointer +char input[ARB] #I Input image name +char output[maxchar] #O Output image name +int maxchar #I Maximum characters in output name. + +int i, fd +pointer sp, section, lv, pv1, pv2 + +int stropen(), strlen() +bool streq() + +begin + call smark (sp) + call salloc (section, SZ_FNAME, TY_CHAR) + call salloc (lv, IM_MAXDIM, TY_LONG) + call salloc (pv1, IM_MAXDIM, TY_LONG) + call salloc (pv2, IM_MAXDIM, TY_LONG) + + # Get endpoint coordinates in original image. + call amovkl (long(1), Meml[lv], IM_MAXDIM) + call aclrl (Meml[pv1], IM_MAXDIM) + call imaplv (im, Meml[lv], Meml[pv1], 2) + call amovl (IM_LEN(im,1), Meml[lv], IM_NDIM(im)) + call aclrl (Meml[pv2], IM_MAXDIM) + call imaplv (im, Meml[lv], Meml[pv2], 2) + + # Set image section. + fd = stropen (Memc[section], SZ_FNAME, NEW_FILE) + call fprintf (fd, "[") + do i = 1, IM_MAXDIM { + if (Meml[pv1+i-1] != Meml[pv2+i-1]) + call fprintf (fd, "*") + else if (Meml[pv1+i-1] != 0) { + call fprintf (fd, "%d") + call pargi (Meml[pv1+i-1]) + } else + break + call fprintf (fd, ",") + } + call close (fd) + i = strlen (Memc[section]) + Memc[section+i-1] = ']' + + if (streq ("[*,*]", Memc[section])) + Memc[section] = EOS + + # Strip existing image section and add new section. + call imgimage (input, output, maxchar) + call strcat (Memc[section], output, maxchar) + +# if (Memc[section] == EOS) +# call imgimage (input, output, maxchar) +# else +# call strcpy (input, output, maxchar) + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/iegnfr.x b/pkg/images/tv/imexamine/iegnfr.x new file mode 100644 index 00000000..0a8fb30d --- /dev/null +++ b/pkg/images/tv/imexamine/iegnfr.x @@ -0,0 +1,61 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include "imexam.h" + +# IE_GETNFRAMES -- Determine the number of image display frames. If the +# display can be accessed at all we assume there is always at least one +# frame; beyond that presence of a valid WCS is used to test whether we +# are interested in looking at a frame. + +int procedure ie_getnframes (ie) + +pointer ie #I imexamine descriptor + +pointer sp, imname, ds, iw +int server, nframes, status, i + +int clgeti(), strncmp(), imd_wcsver() +pointer imd_mapframe(), iw_open() +errchk imd_mapframe, clgeti + +begin + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + + nframes = clgeti ("nframes") + if (nframes == 0) { + # Try to automatically determine the number of frames. + ds = IE_DS(ie) + if (ds == NULL) + ds = imd_mapframe (1, READ_WRITE, NO) + + # If we are talking to a simple image display we assume the device + # has 4 frames (until more general display interfaces come along). + # Servers are more complicated because the number of frames is + # dynamically configurable, even while imexamine is running. + # We use the WCS query to try to count the current number of + # allocated frames in the case of a server device. + + server = IM_LEN(ds,4) + if (server == YES && imd_wcsver() != 0) { + nframes = 1 + do i = 1, MAX_FRAMES { + iferr (iw = iw_open (ds, i, Memc[imname], SZ_FNAME, status)) + next + call iw_close (iw) + if (strncmp (Memc[imname], "[NOSUCHFRAME]", 3) != 0) + nframes = max (nframes, i) + } + } else + nframes = 4 + + if (IE_DS(ie) == NULL) + call imunmap (ds) + } + + IE_NFRAMES(ie) = max (nframes, IE_DFRAME(ie)/100) + call sfree (sp) + + return (nframes) +end diff --git a/pkg/images/tv/imexamine/iegraph.x b/pkg/images/tv/imexamine/iegraph.x new file mode 100644 index 00000000..edfa28c2 --- /dev/null +++ b/pkg/images/tv/imexamine/iegraph.x @@ -0,0 +1,145 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <gset.h> +include "imexam.h" + +define MTYPES "|point|box|plus|cross|circle|hebar|vebar|hline|vline|diamond|" +define IE_GBUF 0.10 # Buffer around data +define IE_SZTITLE 512 # Size of multiline title + + +# IE_GRAPH -- Make a graph +# This procedure is used by most of the different graph types to provide +# consistency in features and parameters. The parameters are read using +# the pset pointer. + +procedure ie_graph (gp, mode, pp, param, x, y, npts, label, format) + +pointer gp # GIO pointer +int mode # Mode +pointer pp # PSET pointer +char param[ARB] # Parameter string +real x[npts] # X data +real y[npts] # Y data +int npts # Number of points +char label # Default x label +char format # Default x format + +int i, marks[10], linepattern, patterns[4], clgpseti(), btoi(), strdic() +pointer sp, title, xlabel, ylabel +real x1, x2, y1, y2, wx1, wx2, wy1, wy2, temp, szmarker +real clgpsetr(), ie_iformatr() +bool clgpsetb(), streq() + +data patterns/GL_SOLID, GL_DASHED, GL_DOTTED, GL_DOTDASH/ +data marks/GM_POINT, GM_BOX, GM_PLUS, GM_CROSS, GM_CIRCLE, GM_HEBAR, + GM_VEBAR, GM_HLINE, GM_VLINE, GM_DIAMOND/ + +begin + call smark (sp) + call salloc (xlabel, SZ_LINE, TY_CHAR) + + # If a new graph setup all the axes and labeling options and then + # make the graph. + + if (mode == NEW_FILE) { + call gclear (gp) + + linepattern = 0 + + x1 = ie_iformatr (clgpsetr (pp, "x1"), format) + x2 = ie_iformatr (clgpsetr (pp, "x2"), format) + y1 = clgpsetr (pp, "y1") + y2 = clgpsetr (pp, "y2") + + if (IS_INDEF (x1) || IS_INDEF (x2)) + call gascale (gp, x, npts, 1) + if (IS_INDEF (y1) || IS_INDEF (y2)) + call gascale (gp, y, npts, 2) + + call gswind (gp, x1, x2, y1, y2) + call ggwind (gp, wx1, wx2, wy1, wy2) + + temp = wx2 - wx1 + if (IS_INDEF (x1)) + wx1 = wx1 - IE_GBUF * temp + if (IS_INDEF (x2)) + wx2 = wx2 + IE_GBUF * temp + + temp = wy2 - wy1 + if (IS_INDEF (y1)) + wy1 = wy1 - IE_GBUF * temp + if (IS_INDEF (y2)) + wy2 = wy2 + IE_GBUF * temp + + call gswind (gp, wx1, wx2, wy1, wy2) + call gsetr (gp, G_ASPECT, 0.) + call gseti (gp, G_ROUND, btoi (clgpsetb (pp, "round"))) + + i = GW_LINEAR + if (clgpsetb (pp, "logx")) + i = GW_LOG + call gseti (gp, G_XTRAN, i) + i = GW_LINEAR + if (clgpsetb (pp, "logy")) + i = GW_LOG + call gseti (gp, G_YTRAN, i) + + if (clgpsetb (pp, "box")) { + # Get number of major and minor tick marks. + call gseti (gp, G_XNMAJOR, clgpseti (pp, "majrx")) + call gseti (gp, G_XNMINOR, clgpseti (pp, "minrx")) + call gseti (gp, G_YNMAJOR, clgpseti (pp, "majry")) + call gseti (gp, G_YNMINOR, clgpseti (pp, "minry")) + + # Label tick marks on axes? + call gsets (gp, G_XTICKFORMAT, format) + call gseti (gp, G_LABELTICKS, + btoi (clgpsetb (pp, "ticklabels"))) + + # Fetch labels and plot title string. + call salloc (title, IE_SZTITLE, TY_CHAR) + call salloc (ylabel, SZ_LINE, TY_CHAR) + + if (clgpsetb (pp, "banner")) { + call sysid (Memc[title], IE_SZTITLE) + call strcat ("\n", Memc[title], IE_SZTITLE) + call strcat (param, Memc[title], IE_SZTITLE) + } else + Memc[title] = EOS + + call clgpset (pp, "title", Memc[xlabel], SZ_LINE) + if (Memc[xlabel] != EOS) { + call strcat ("\n", Memc[title], IE_SZTITLE) + call strcat (Memc[xlabel], Memc[title], IE_SZTITLE) + } + call clgpset (pp, "xlabel", Memc[xlabel], SZ_LINE) + call clgpset (pp, "ylabel", Memc[ylabel], SZ_LINE) + + if (streq ("wcslabel", Memc[xlabel])) + call strcpy (label, Memc[xlabel], SZ_LINE) + + call glabax (gp, Memc[title], Memc[xlabel], Memc[ylabel]) + } + } + + # Draw the data. + if (clgpsetb (pp, "pointmode")) { + call clgpset (pp, "marker", Memc[xlabel], SZ_LINE) + i = strdic (Memc[xlabel], Memc[xlabel], SZ_LINE, MTYPES) + if (i == 0) + i = 2 + if (marks[i] == GM_POINT) + szmarker = 0.0 + else + szmarker = clgpsetr (pp, "szmarker") + call gpmark (gp, x, y, npts, marks[i], szmarker, szmarker) + } else { + linepattern = min (4, linepattern + 1) + call gseti (gp, G_PLTYPE, patterns[linepattern]) + call gpline (gp, x, y, npts) + } + call gflush (gp) + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/iehimexam.x b/pkg/images/tv/imexamine/iehimexam.x new file mode 100644 index 00000000..4a0fd150 --- /dev/null +++ b/pkg/images/tv/imexamine/iehimexam.x @@ -0,0 +1,193 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include "imexam.h" + +define HGM_TYPES "|line|box|" +define HGM_LINE 1 # line vectors for histogram plot +define HGM_BOX 2 # box vectors for histogram plot + +# IE_HIMEXAM -- Compute and plot or list a histogram. +# If the GIO pointer is NULL list the histogram otherwise make a graph. + +procedure ie_himexam (gp, mode, ie, x, y) + +pointer gp # GIO pointer (NULL for histogram listing) +int mode # Mode +pointer ie # Structure pointer +real x, y # Center coordinate + +real z1, z2, dz, zmin, zmax +int i, j, x1, x2, y1, y2, nx, ny, npts, nbins, nbins1, nlevels, nwide +pointer pp, sp, hgm, title, im, data, xp, yp + +int clgpseti() +real clgpsetr() +bool clgpsetb(), fp_equalr() +pointer clopset(), ie_gimage(), ie_gdata() + +begin + # Get the image and return on error. + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + # Use last graph coordinate if redrawing. Close last graph pset + # pointer if making new graph. + + if (gp != NULL) { + if (!IS_INDEF(x)) + IE_X1(ie) = x + if (!IS_INDEF(y)) + IE_Y1(ie) = y + + z1 = IE_X1(ie) + z2 = IE_Y1(ie) + + if (IE_PP(ie) != NULL) + call clcpset (IE_PP(ie)) + } else { + z1 = x + z2 = y + } + + # Get the data. + pp = clopset ("himexam") + nx = clgpseti (pp, "ncolumns") + ny = clgpseti (pp, "nlines") + x1 = z1 - (nx - 1) / 2 + 0.5 + x2 = z1 + nx / 2 + 0.5 + y1 = z2 - (ny - 1) / 2 + 0.5 + y2 = z2 + ny / 2 + 0.5 + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + return + } + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + npts = nx * ny + + # Get default histogram resolution. + nbins = clgpseti (pp, "nbins") + + # Get histogram range. + z1 = clgpsetr (pp, "z1") + z2 = clgpsetr (pp, "z2") + + # Use data limits for INDEF limits. + if (IS_INDEFR(z1) || IS_INDEFR(z2)) { + call alimr (Memr[data], npts, zmin, zmax) + if (IS_INDEFR(z1)) + z1 = zmin + if (IS_INDEFR(z2)) + z2 = zmax + } + + if (z1 > z2) { + dz = z1; z1 = z2; z2 = dz + } + + # Adjust the resolution of the histogram and/or the data range + # so that an integral number of data values map into each + # histogram bin (to avoid aliasing effects). + + if (clgpsetb (pp, "autoscale")) { + switch (IM_PIXTYPE(im)) { + case TY_SHORT, TY_USHORT, TY_INT, TY_LONG: + nlevels = nint (z2) - nint (z1) + nwide = max (1, nint (real (nlevels) / real (nbins))) + nbins = max (1, nint (real (nlevels) / real (nwide))) + z2 = nint (z1) + nbins * nwide + } + } + + # Test for constant valued image, which causes zero divide in ahgm. + if (fp_equalr (z1, z2)) { + call eprintf ("Warning: Image `%s' has no data range.\n") + call pargstr (IE_IMAGE(ie)) + return + } + + # The extra bin counts the pixels that equal z2 and shifts the + # remaining bins to evenly cover the interval [z1,z2]. + # Note that real numbers could be handled better - perhaps + # adjust z2 upward by ~ EPSILONR (in ahgm itself). + + nbins1 = nbins + 1 + + # Initialize the histogram buffer and image line vector. + call smark (sp) + call salloc (hgm, nbins1, TY_INT) + call aclri (Memi[hgm], nbins1) + + call ahgmr (Memr[data], npts, Memi[hgm], nbins1, z1, z2) + + # "Correct" the topmost bin for pixels that equal z2. Each + # histogram bin really wants to be half open. + + if (clgpsetb (pp, "top_closed")) + Memi[hgm+nbins-1] = Memi[hgm+nbins-1] + Memi[hgm+nbins1-1] + + # List or plot the histogram. In list format, the bin value is the + # z value of the left side (start) of the bin. + + dz = (z2 - z1) / real (nbins) + + if (gp != NULL) { + # Draw the plot. + if (clgpsetb (pp, "pointmode")) { + nbins1 = nbins + call salloc (xp, nbins1, TY_REAL) + call salloc (yp, nbins1, TY_REAL) + call achtir (Memi[hgm], Memr[yp], nbins1) + Memr[xp] = z1 + dz / 2. + do i = 1, nbins1 - 1 + Memr[xp+i] = Memr[xp+i-1] + dz + } else { + nbins1 = 2 * nbins + call salloc (xp, nbins1, TY_REAL) + call salloc (yp, nbins1, TY_REAL) + Memr[xp] = z1 + Memr[yp] = Memi[hgm] + j = 0 + do i = 1, nbins - 1 { + Memr[xp+j+1] = Memr[xp+j] + dz + Memr[yp+j+1] = Memr[yp+j] + j = j + 1 + Memr[xp+j+1] = Memr[xp+j] + Memr[yp+j+1] = Memi[hgm+i] + j = j + 1 + } + Memr[xp+j+1] = Memr[xp+j] + dz + Memr[yp+j+1] = Memr[yp+j] + } + + call salloc (title, IE_SZTITLE, TY_CHAR) + call sprintf (Memc[title], IE_SZTITLE, + "%s[%d:%d,%d:%d]: Histogram from z1=%g to z2=%g, nbins=%d\n%s") + call pargstr (IE_IMNAME(ie)) + call pargi (x1) + call pargi (x2) + call pargi (y1) + call pargi (y2) + call pargr (z1) + call pargr (z2) + call pargi (nbins) + call pargstr (IM_TITLE(im)) + call ie_graph (gp, mode, pp, Memc[title], Memr[xp], + Memr[yp], nbins1, "", "") + + IE_PP(ie) = pp + } else { + do i = 1, nbins { + call printf ("%g %d\n") + call pargr (z1 + (i-1) * dz) + call pargi (Memi[hgm+i-1]) + } + call clcpset (pp) + } + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/ieimname.x b/pkg/images/tv/imexamine/ieimname.x new file mode 100644 index 00000000..3b1bd5e9 --- /dev/null +++ b/pkg/images/tv/imexamine/ieimname.x @@ -0,0 +1,33 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +# IE_IMNAME -- Get the name of the image displayed in a display frame. + +procedure ie_imname (ds, frame, imname, maxch) + +pointer ds #I display descriptor +int frame #I display frame +char imname[maxch] #O image name +int maxch #I max chars out + +int snx, sny, dx, dy, dnx, dny, status, imd_query_map() +real sx, sy +pointer sp, reg, dname, iw +pointer iw_open() +errchk imd_query_map, iw_open + +begin + call smark (sp) + call salloc (reg, SZ_FNAME, TY_CHAR) + call salloc (dname, SZ_FNAME, TY_CHAR) + + if (imd_query_map (frame, Memc[reg], sx, sy, snx, sny, dx, dy, dnx, dny, + Memc[dname]) == ERR) { + iw = iw_open (ds, frame/100, Memc[dname], SZ_FNAME, status) + call iw_close (iw) + } + + # call imgimage (Memc[dname], imname, maxch) + call strcpy (Memc[dname], imname, maxch) + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/iejimexam.x b/pkg/images/tv/imexamine/iejimexam.x new file mode 100644 index 00000000..46a4c910 --- /dev/null +++ b/pkg/images/tv/imexamine/iejimexam.x @@ -0,0 +1,473 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include <gset.h> +include <mach.h> +include "imexam.h" + + +# IE_JIMEXAM -- 1D profile plot and gaussian fit parameters. +# If no GIO pointer is given then only the fit parameters are printed. +# The fitting uses a Levenberg-Marquardt nonlinear chi square minimization. + +procedure ie_jimexam (gp, mode, ie, x, y, axis) + +pointer gp +pointer ie +int mode +real x, y +int axis + +int navg, order, clgpseti() +bool center, background, clgpsetb() +real sigma, width, rplot, clgpsetr() + +int i, j, k, nx, ny, x1, x2, y1, y2, nfit, flag[5] +real xc, yc, bkg, r, dr, fit[5], xfit, yfit, asumr(), amedr() +pointer sp, title, avstr, im, pp, data, xs, ys, ptr +pointer clopset(), ie_gimage(), ie_gdata() + +errchk ie_gdata, mr_solve + +begin + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + # Get parameters + if (IE_PP(ie) != NULL) + call clcpset (IE_PP(ie)) + if (axis == 1) + IE_PP(ie) = clopset ("jimexam") + else + IE_PP(ie) = clopset ("kimexam") + pp = IE_PP(ie) + navg = clgpseti (pp, "naverage") + center = clgpsetb (pp, "center") + background = clgpsetb (pp, "background") + sigma = clgpsetr (pp, "sigma") + rplot = clgpsetr (pp, "rplot") + if (background) { + order = clgpsetr (pp, "xorder") + width = clgpsetr (pp, "width") + } + + # If the initial center is INDEF then use the previous value. + if (!IS_INDEF(x)) + IE_X1(ie) = x + if (!IS_INDEF(y)) + IE_Y1(ie) = y + + if (axis == 1) { + xc = IE_X1(ie) + yc = IE_Y1(ie) + } else { + xc = IE_Y1(ie) + yc = IE_X1(ie) + } + + # Get data + r = max (rplot, 8 * sigma + width) + x1 = xc - r + x2 = xc + r + y1 = nint (yc) - (navg - 1) / 2 + y2 = nint (yc) + navg / 2 + iferr { + if (axis == 1) + data = ie_gdata (im, x1, x2, y1, y2) + else + data = ie_gdata (im, y1, y2, x1, x2) + } then { + call erract (EA_WARN) + return + } + + # Compute average vector + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + yc = (y1 + y2) / 2. + + call smark (sp) + call salloc (xs, nx, TY_REAL) + call salloc (ys, nx, TY_REAL) + call salloc (title, IE_SZTITLE, TY_CHAR) + call salloc (avstr, SZ_LINE, TY_CHAR) + + ptr = data + if (axis == 1) { + call sprintf (Memc[avstr], SZ_LINE, "Lines %d-%d") + call pargi (y1) + call pargi (y2) + call amovr (Memr[ptr], Memr[ys], nx) + ptr = ptr + nx + do i = 2, ny { + call aaddr (Memr[ptr], Memr[ys], Memr[ys], nx) + ptr = ptr + nx + } + call adivkr (Memr[ys], real (ny), Memr[ys], nx) + } else { + call sprintf (Memc[avstr], SZ_LINE, "Columns %d-%d") + call pargi (y1) + call pargi (y2) + do i = 0, nx-1 { + Memr[ys+i] = asumr (Memr[ptr], ny) / ny + ptr = ptr + ny + } + } + + # Set default background + bkg = 0. + if (background) { + r = 4 * sigma + ptr = xs + do i = 0, nx-1 { + if (abs (xc - x1 - i) > r) { + Memr[ptr] = Memr[ys+i] + ptr = ptr + 1 + } + } + if (ptr > xs) + bkg = amedr (Memr[xs], ptr-xs) + } + + # Convert to WCS + if (axis == 1) { + call ie_mwctran (ie, xc, yc, xfit, yfit) + call ie_mwctran (ie, xc+sigma, yc, r, yfit) + dr = abs (xfit - r) + do i = 0, nx-1 + call ie_mwctran (ie, real(x1+i), yc, Memr[xs+i], yfit) + } else { + call ie_mwctran (ie, yc, xc, yfit, xfit) + call ie_mwctran (ie, yc, xc+sigma, yfit, r) + dr = abs (xfit - r) + do i = 0, nx-1 + call ie_mwctran (ie, yc, real(x1+i), yfit, Memr[xs+i]) + } + + # Set initial fit parameters + k = max (0, nint (xc - x1)) + fit[1] = bkg + fit[2] = 0. + fit[3] = Memr[ys+k] - fit[1] + fit[4] = xfit + fit[5] = dr + + # Do fitting. + nfit = 1 + flag[1] = 3 + + # Add centering if desired + if (center) { + nfit = nfit + 1 + flag[nfit] = 4 + call ie_gfit (Memr[xs], Memr[ys], nx, fit, flag, nfit) + } + + # Add sigma + nfit = nfit + 1 + flag[nfit] = 5 + call ie_gfit (Memr[xs], Memr[ys], nx, fit, flag, nfit) + + # Now add background if desired + if (background) { + if (order == 1) { + nfit = nfit + 1 + flag[nfit] = 1 + call ie_gfit (Memr[xs], Memr[ys], nx, fit, flag, nfit) + } else if (order == 2) { + nfit = nfit + 2 + flag[nfit-1] = 1 + flag[nfit] = 2 + call ie_gfit (Memr[xs], Memr[ys], nx, fit, flag, nfit) + } + } + + # Plot the profile and overplot the gaussian fit. + call sprintf (Memc[title], IE_SZTITLE, "%s: %s\n%s") + call pargstr (IE_IMNAME(ie)) + call pargstr (Memc[avstr]) + call pargstr (IM_TITLE(im)) + + j = max (0, int (xc - x1 - rplot)) + k = min (nx-1, nint (xc - x1 + rplot)) + if (axis == 1) + call ie_graph (gp, mode, pp, Memc[title], + Memr[xs+j], Memr[ys+j], k-j+1, IE_XLABEL(ie), IE_XFORMAT(ie)) + else + call ie_graph (gp, mode, pp, Memc[title], + Memr[xs+j], Memr[ys+j], k-j+1, IE_YLABEL(ie), IE_YFORMAT(ie)) + + call gseti (gp, G_PLTYPE, 2) + xfit = min (Memr[xs+j], Memr[xs+k]) + r = (xfit - fit[4]) / fit[5] + dr = abs ((Memr[xs+k] - Memr[xs+j]) / (k - j)) + if (abs (r) < 7.) + yfit = fit[1] + fit[2] * xfit + fit[3] * exp (-r**2 / 2.) + else + yfit = fit[1] + fit[2] * xfit + call gamove (gp, xfit, yfit) + repeat { + xfit = xfit + 0.2 * dr + r = (xfit - fit[4]) / fit[5] + if (abs (r) < 7.) + yfit = fit[1] + fit[2] * xfit + fit[3] * exp (-r**2 / 2.) + else + yfit = fit[1] + fit[2] * xfit + call gadraw (gp, xfit, yfit) + } until (xfit >= max (Memr[xs+j], Memr[xs+k])) + call gseti (gp, G_PLTYPE, 1) + + # Print the fit values + call printf ("%s: center=%7g peak=%7g sigma=%7.4g fwhm=%7.4g bkg=%7g\n") + call pargstr (Memc[avstr]) + call pargr (fit[4]) + call pargr (fit[3]) + call pargr (fit[5]) + call pargr (2.35482*fit[5]) + call pargr (fit[1]+fit[2]*fit[4]) + + if (IE_LOGFD(ie) != NULL) { + call fprintf (IE_LOGFD(ie), + "%s: center=%7g peak=%7g sigma=%5.3f fwhm=%5.3f bkg=%7g\n") + call pargstr (Memc[avstr]) + call pargr (fit[4]) + call pargr (fit[3]) + call pargr (fit[5]) + call pargr (2.35482*fit[5]) + call pargr (fit[1]+fit[2]*fit[4]) + } + + call sfree (sp) +end + + +# IE_GFIT -- 1D Gaussian fit. + +procedure ie_gfit (xs, ys, nx, fit, flag, nfit) + +real xs[nx], ys[nx] # Vector to be fit +int nx # Number of points +real fit[5] # Fit parameters +int flag[nfit] # Flag for parameters to be fit +int nfit # Number of parameters to be fit + +int i +real chi1, chi2, mr + +begin + chi2 = MAX_REAL + mr = -1. + i = 0 + repeat { + call mr_solve (xs, ys, nx, fit, flag, 5, nfit, mr, chi1) + if (chi2 - chi1 > 1.) + i = 0 + else + i = i + 1 + chi2 = chi1 + } until (i == 3) + mr = 0. + call mr_solve (xs, ys, nx, fit, flag, 5, nfit, mr, chi1) + + fit[5] = abs (fit[5]) +end + + +# DERIVS -- Compute model and derivatives for MR_SOLVE procedure. +# +# I(x) = A1 + A2 * x + A3 exp {-[(x - A4) / A5]**2 / 2.} +# +# where the params are A1-A5. + +procedure derivs (x, a, y, dyda, na) + +real x # X value to be evaluated +real a[na] # Parameters +real y # Function value +real dyda[na] # Derivatives +int na # Number of parameters + +real arg, ex, fac + +begin + arg = (x - a[4]) / a[5] + if (abs (arg) < 7.) + ex = exp (-arg**2 / 2.) + else + ex = 0. + fac = a[3] * ex * arg + + y = a[1] + a[2] * x + a[3] * ex + + dyda[1] = 1. + dyda[2] = x + dyda[3] = ex + dyda[4] = fac / a[5] + dyda[5] = fac * arg / a[5] +end + + +# MR_SOLVE -- Levenberg-Marquardt nonlinear chi square minimization. +# +# Use the Levenberg-Marquardt method to minimize the chi squared of a set +# of paraemters. The parameters being fit are indexed by the flag array. +# To initialize the Marquardt parameter, MR, is less than zero. After that +# the parameter is adjusted as needed. To finish set the parameter to zero +# to free memory. This procedure requires a subroutine, DERIVS, which +# takes the derivatives of the function being fit with respect to the +# parameters. There is no limitation on the number of parameters or +# data points. For a description of the method see NUMERICAL RECIPES +# by Press, Flannery, Teukolsky, and Vetterling, p523. + +procedure mr_solve (x, y, npts, params, flags, np, nfit, mr, chisq) + +real x[npts] # X data array +real y[npts] # Y data array +int npts # Number of data points +real params[np] # Parameter array +int flags[np] # Flag array indexing parameters to fit +int np # Number of parameters +int nfit # Number of parameters to fit +real mr # MR parameter +real chisq # Chi square of fit + +int i +real chisq1 +pointer new, a1, a2, delta1, delta2 + +errchk mr_invert + +begin + # Allocate memory and initialize. + if (mr < 0.) { + call mfree (new, TY_REAL) + call mfree (a1, TY_REAL) + call mfree (a2, TY_REAL) + call mfree (delta1, TY_REAL) + call mfree (delta2, TY_REAL) + + call malloc (new, np, TY_REAL) + call malloc (a1, nfit*nfit, TY_REAL) + call malloc (a2, nfit*nfit, TY_REAL) + call malloc (delta1, nfit, TY_REAL) + call malloc (delta2, nfit, TY_REAL) + + call amovr (params, Memr[new], np) + call mr_eval (x, y, npts, Memr[new], flags, np, Memr[a2], + Memr[delta2], nfit, chisq) + mr = 0.001 + } + + # Restore last good fit and apply the Marquardt parameter. + call amovr (Memr[a2], Memr[a1], nfit * nfit) + call amovr (Memr[delta2], Memr[delta1], nfit) + do i = 1, nfit + Memr[a1+(i-1)*(nfit+1)] = Memr[a2+(i-1)*(nfit+1)] * (1. + mr) + + # Matrix solution. + call mr_invert (Memr[a1], Memr[delta1], nfit) + + # Compute the new values and curvature matrix. + do i = 1, nfit + Memr[new+flags[i]-1] = params[flags[i]] + Memr[delta1+i-1] + call mr_eval (x, y, npts, Memr[new], flags, np, Memr[a1], + Memr[delta1], nfit, chisq1) + + # Check if chisq has improved. + if (chisq1 < chisq) { + mr = max (EPSILONR, 0.1 * mr) + chisq = chisq1 + call amovr (Memr[a1], Memr[a2], nfit * nfit) + call amovr (Memr[delta1], Memr[delta2], nfit) + call amovr (Memr[new], params, np) + } else + mr = 10. * mr + + if (mr == 0.) { + call mfree (new, TY_REAL) + call mfree (a1, TY_REAL) + call mfree (a2, TY_REAL) + call mfree (delta1, TY_REAL) + call mfree (delta2, TY_REAL) + } +end + + +# MR_EVAL -- Evaluate curvature matrix. This calls procedure DERIVS. + +procedure mr_eval (x, y, npts, params, flags, np, a, delta, nfit, chisq) + +real x[npts] # X data array +real y[npts] # Y data array +int npts # Number of data points +real params[np] # Parameter array +int flags[np] # Flag array indexing parameters to fit +int np # Number of parameters +real a[nfit,nfit] # Curvature matrix +real delta[nfit] # Delta array +int nfit # Number of parameters to fit +real chisq # Chi square of fit + +int i, j, k +real ymod, dy, dydpj, dydpk +pointer sp, dydp + +begin + call smark (sp) + call salloc (dydp, np, TY_REAL) + + do j = 1, nfit { + do k = 1, j + a[j,k] = 0. + delta[j] = 0. + } + + chisq = 0. + do i = 1, npts { + call derivs (x[i], params, ymod, Memr[dydp], np) + dy = y[i] - ymod + do j = 1, nfit { + dydpj = Memr[dydp+flags[j]-1] + delta[j] = delta[j] + dy * dydpj + do k = 1, j { + dydpk = Memr[dydp+flags[k]-1] + a[j,k] = a[j,k] + dydpj * dydpk + } + } + chisq = chisq + dy * dy + } + + do j = 2, nfit + do k = 1, j-1 + a[k,j] = a[j,k] + + call sfree (sp) +end + + +# MR_INVERT -- Solve a set of linear equations using Householder transforms. + +procedure mr_invert (a, b, n) + +real a[n,n] # Input matrix and returned inverse +real b[n] # Input RHS vector and returned solution +int n # Dimension of input matrices + +int krank +real rnorm +pointer sp, h, g, ip + +begin + call smark (sp) + call salloc (h, n, TY_REAL) + call salloc (g, n, TY_REAL) + call salloc (ip, n, TY_INT) + + call hfti (a, n, n, n, b, n, 1, 1E-10, krank, rnorm, + Memr[h], Memr[g], Memi[ip]) + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/ielimexam.x b/pkg/images/tv/imexamine/ielimexam.x new file mode 100644 index 00000000..9b1c490d --- /dev/null +++ b/pkg/images/tv/imexamine/ielimexam.x @@ -0,0 +1,81 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include "imexam.h" + + +# IE_LIMEXAM -- Make a line plot +# If the line is INDEF then use the last line. + +procedure ie_limexam (gp, mode, ie, y) + +pointer gp # GIO pointer +int mode # Mode +pointer ie # Structure pointer +real y # Line + +real yavg, junk +int i, x1, x2, y1, y2, nx, ny, npts +pointer sp, title, im, data, ptr, xp, yp + +int clgpseti() +pointer clopset(), ie_gimage(), ie_gdata() + +begin + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + if (IE_PP(ie) != NULL) + call clcpset (IE_PP(ie)) + IE_PP(ie) = clopset ("limexam") + + if (!IS_INDEF(y)) + IE_Y1(ie) = y + + ny = clgpseti (IE_PP(ie), "naverage") + x1 = INDEFI + x2 = INDEFI + y1 = IE_Y1(ie) - (ny - 1) / 2 + 0.5 + y2 = IE_Y1(ie) + ny / 2 + 0.5 + yavg = (y1 + y2) / 2. + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + return + } + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + npts = nx * ny + + call smark (sp) + call salloc (title, IE_SZTITLE, TY_CHAR) + call salloc (xp, nx, TY_REAL) + + do i = 1, nx + call ie_mwctran (ie, real(i), yavg, Memr[xp+i-1], junk) + + if (ny > 1) { + ptr = data + call salloc (yp, nx, TY_REAL) + call amovr (Memr[ptr], Memr[yp], nx) + do i = 2, ny { + ptr = ptr + nx + call aaddr (Memr[ptr], Memr[yp], Memr[yp], nx) + } + call adivkr (Memr[yp], real (ny), Memr[yp], nx) + } else + yp = data + + call sprintf (Memc[title], IE_SZTITLE, "%s: Lines %d - %d\n%s") + call pargstr (IE_IMNAME(ie)) + call pargi (y1) + call pargi (y2) + call pargstr (IM_TITLE(im)) + + call ie_graph (gp, mode, IE_PP(ie), Memc[title], Memr[xp], + Memr[yp], nx, IE_XLABEL(ie), IE_XFORMAT(ie)) + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/iemw.x b/pkg/images/tv/imexamine/iemw.x new file mode 100644 index 00000000..185cfbaa --- /dev/null +++ b/pkg/images/tv/imexamine/iemw.x @@ -0,0 +1,191 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <mwset.h> +include "imexam.h" + + +# IE_MWINIT -- Initialize MWCS + +procedure ie_mwinit (ie) + +pointer ie # IMEXAM descriptor + +int i, j, wcsdim, mw_stati(), nowhite(), stridxs() +pointer im, mw, ctlw, ctwl, mw_openim(), mw_sctran() +pointer sp, axno, axval, str1, str2 +bool streq() +errchk mw_openim, mw_sctran + +begin + im = IE_IM(ie) + mw = IE_MW(ie) + + if (mw != NULL) { + call mw_close (mw) + IE_MW(ie) = mw + } + + IE_XLABEL(ie) = EOS + IE_YLABEL(ie) = EOS + call clgstr ("xformat", IE_XFORMAT(ie), IE_SZFORMAT) + call clgstr ("yformat", IE_YFORMAT(ie), IE_SZFORMAT) + i = nowhite (IE_XFORMAT(ie), IE_XFORMAT(ie), IE_SZFORMAT) + i = nowhite (IE_YFORMAT(ie), IE_YFORMAT(ie), IE_SZFORMAT) + + if (im == NULL || im == IE_DS(ie)) + return + + call smark (sp) + call salloc (axno, IM_MAXDIM, TY_INT) + call salloc (axval, IM_MAXDIM, TY_INT) + call salloc (str1, SZ_LINE, TY_CHAR) + call salloc (str2, SZ_LINE, TY_CHAR) + + mw = mw_openim (im) + call mw_seti (mw, MW_USEAXMAP, NO) + wcsdim = mw_stati (mw, MW_NDIM) + call mw_gaxmap (mw, Memi[axno], Memi[axval], wcsdim) + IE_P1(ie) = 1 + IE_P2(ie) = 2 + do i = 1, wcsdim { + j = Memi[axno+i-1] + if (j == 0) + IE_IN(ie,i) = 1 + else if (j == 1) + IE_P1(ie) = i + else if (j == 2) + IE_P2(ie) = i + } + ctlw = mw_sctran (mw, "logical", IE_WCSNAME(ie), 0) + ctwl = mw_sctran (mw, IE_WCSNAME(ie), "logical", 0) + + # Set coordinate labels and formats + i = IE_P1(ie) + j = IE_P2(ie) + if (streq (IE_WCSNAME(ie), "logical")) { + call strcpy ("Column (pixels)", IE_XLABEL(ie), IE_SZFNAME) + call strcpy ("Line (pixels)", IE_YLABEL(ie), IE_SZFNAME) + } else if (streq (IE_WCSNAME(ie), "physical")) { + if (i == 1) + call strcpy ("Column (pixels)", IE_XLABEL(ie), IE_SZFNAME) + else if (i == 2) + call strcpy ("Line (pixels)", IE_XLABEL(ie), IE_SZFNAME) + else + call strcpy ("Pixels", IE_XLABEL(ie), IE_SZFNAME) + if (j == 1) + call strcpy ("Column (pixels)", IE_YLABEL(ie), IE_SZFNAME) + else if (j == 2) + call strcpy ("Line (pixels)", IE_YLABEL(ie), IE_SZFNAME) + else + call strcpy ("Pixels", IE_YLABEL(ie), IE_SZFNAME) + } else { + ifnoerr (call mw_gwattrs (mw, i, "label", Memc[str1], SZ_LINE)) { + ifnoerr (call mw_gwattrs (mw, i, "units", Memc[str2],SZ_LINE)) { + call sprintf (IE_XLABEL(ie), IE_SZFNAME, "%s (%s)") + call pargstr (Memc[str1]) + call pargstr (Memc[str2]) + } else { + call sprintf (IE_XLABEL(ie), IE_SZFNAME, "%s") + call pargstr (Memc[str1]) + } + } + if (IE_XFORMAT(ie) != '%') + ifnoerr (call mw_gwattrs (mw, i, "format", Memc[str1], SZ_LINE)) + call strcpy (Memc[str1], IE_XFORMAT(ie), IE_SZFORMAT) + + ifnoerr (call mw_gwattrs (mw, j, "label", Memc[str1], SZ_LINE)) { + ifnoerr (call mw_gwattrs (mw, j, "units", Memc[str2],SZ_LINE)) { + call sprintf (IE_YLABEL(ie), IE_SZFNAME, "%s (%s)") + call pargstr (Memc[str1]) + call pargstr (Memc[str2]) + } else { + call sprintf (IE_YLABEL(ie), IE_SZFNAME, "%s") + call pargstr (Memc[str1]) + } + } + if (IE_YFORMAT(ie) != '%') + ifnoerr (call mw_gwattrs (mw, j, "format", Memc[str1], SZ_LINE)) + call strcpy (Memc[str1], IE_YFORMAT(ie), IE_SZFORMAT) + + # Check for equitorial coordinate and reversed formats. + ifnoerr (call mw_gwattrs (mw, i, "axtype", Memc[str1], SZ_LINE)) + if ((streq(Memc[str1],"ra")&&stridxs("hm",IE_XFORMAT(ie))>0) || + (streq(Memc[str1],"dec")&&stridxs("HM",IE_XFORMAT(ie))>0)) { + call strcpy (IE_XFORMAT(ie), Memc[str1], IE_SZFORMAT) + call strcpy (IE_YFORMAT(ie), IE_XFORMAT(ie),IE_SZFORMAT) + call strcpy (Memc[str1], IE_YFORMAT(ie), IE_SZFORMAT) + } + } + + IE_MW(ie) = mw + IE_CTLW(ie) = ctlw + IE_CTWL(ie) = ctwl + IE_WCSDIM(ie) = wcsdim + + call sfree (sp) +end + + +# IE_MWCTRAN -- Evaluate MWCS coordinate + +procedure ie_mwctran (ie, xin, yin, xout, yout) + +pointer ie # IMEXAM descriptor +real xin, yin # Input coordinate +real xout, yout # Output coordinate + +begin + if (IE_MW(ie) == NULL) { + xout = xin + yout = yin + return + } + + IE_IN(ie,IE_P1(ie)) = xin + IE_IN(ie,IE_P2(ie)) = yin + call mw_ctranr (IE_CTLW(ie), IE_IN(ie,1), IE_OUT(ie,1), IE_WCSDIM(ie)) + xout = IE_OUT(ie,IE_P1(ie)) + yout = IE_OUT(ie,IE_P2(ie)) +end + + +# IE_IMWCTRAN -- Evaluate inverse MWCS coordinate + +procedure ie_imwctran (ie, xin, yin, xout, yout) + +pointer ie # IMEXAM descriptor +real xin, yin # Input coordinate +real xout, yout # Output coordinate + +begin + if (IE_MW(ie) == NULL) { + xout = xin + yout = yin + return + } + + IE_OUT(ie,IE_P1(ie)) = xin + IE_OUT(ie,IE_P2(ie)) = yin + call mw_ctranr (IE_CTWL(ie), IE_OUT(ie,1), IE_IN(ie,1), IE_WCSDIM(ie)) + xout = IE_IN(ie,IE_P1(ie)) + yout = IE_IN(ie,IE_P2(ie)) +end + + +# IE_IFORMATR -- Determine the inverse formatted real value +# This temporary routine is used to account for scaling of the H and M formats. + +real procedure ie_iformatr (value, format) + +real value # Value to be inverse formated +char format[ARB] # Format + +int strldxs() + +begin + if (!IS_INDEF(value) && strldxs ("HM", format) > 0) + return (value * 15.) + else + return (value) +end diff --git a/pkg/images/tv/imexamine/ieopenlog.x b/pkg/images/tv/imexamine/ieopenlog.x new file mode 100644 index 00000000..08f754f9 --- /dev/null +++ b/pkg/images/tv/imexamine/ieopenlog.x @@ -0,0 +1,39 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include "imexam.h" + + +# IE_OPENLOG -- Open the log file. + +procedure ie_openlog (ie) + +pointer ie #I imexamine descriptor + +int nowhite(), open() +errchk open, close + +begin + if (IE_LOGFD(ie) != NULL) { + call close (IE_LOGFD(ie)) + IE_LOGFD(ie) = NULL + } + + if (nowhite (IE_LOGFILE(ie), IE_LOGFILE(ie), SZ_FNAME) > 0) { + iferr { + IE_LOGFD(ie) = open (IE_LOGFILE(ie), APPEND, TEXT_FILE) + call printf ("Log file %s open\n") + call pargstr (IE_LOGFILE(ie)) + + if (IE_IM(ie) != NULL) { + call fprintf (IE_LOGFD(ie), "# [%d] %s - %s\n") + call pargi (IE_INDEX(ie)) + call pargstr (IE_IMNAME(ie)) + call pargstr (IM_TITLE(IE_IM(ie))) + } + + } then + call erract (EA_WARN) + } +end diff --git a/pkg/images/tv/imexamine/iepos.x b/pkg/images/tv/imexamine/iepos.x new file mode 100644 index 00000000..7253816b --- /dev/null +++ b/pkg/images/tv/imexamine/iepos.x @@ -0,0 +1,180 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <math.h> +include "imexam.h" + +# IE_POS -- Print cursor position and pixel value or set new origin. +# If the origin is not (0,0) print additional fields. + +procedure ie_pos (ie, x, y, key) + +pointer ie # IMEXAM structure +real x, y # Center of box +int key # Key ('x' positions, 'y' origin) + +pointer im, data +real dx, dy, r, t, wx, wy, xo, yo +int x1, x2, y1, y2 +pointer ie_gimage(), ie_gdata() + +begin + switch (key) { + case 'x': # Print position and pixel value + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + x1 = x + 0.5 + x2 = x + 0.5 + y1 = y + 0.5 + y2 = y + 0.5 + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + return + } + + call printf ("%7.2f %7.2f %7g") + call pargr (x) + call pargr (y) + call pargr (Memr[data]) + + # Print additional fields + if (IE_XORIGIN(ie) != 0. || IE_YORIGIN(ie) != 0.) { + dx = x - IE_XORIGIN(ie) + dy = y - IE_YORIGIN(ie) + r = sqrt (dx * dx + dy * dy) + t = mod (360. + RADTODEG (atan2 (dy, dx)), 360.) + call printf (" %7.f %7.2f %7.2f %7.2f %7.2f %5.1f") + call pargr (IE_XORIGIN(ie)) + call pargr (IE_YORIGIN(ie)) + call pargr (dx) + call pargr (dy) + call pargr (r) + call pargr (t) + } + call printf ("\n") + case 'y': # Set new origin + IE_XORIGIN(ie) = x + IE_YORIGIN(ie) = y + call printf ("Origin: %.2f %.2f\n") + call pargr (IE_XORIGIN(ie)) + call pargr (IE_YORIGIN(ie)) + } + + # Print to logfile if needed. + if (IE_LOGFD(ie) != NULL) { + switch (key) { + case 'x': + call fprintf (IE_LOGFD(ie), "%7.2f %7.2f %7g") + call pargr (x) + call pargr (y) + call pargr (Memr[data]) + if (IE_XORIGIN(ie) != 0. || IE_YORIGIN(ie) != 0.) { + dx = x - IE_XORIGIN(ie) + dy = y - IE_YORIGIN(ie) + r = sqrt (dx * dx + dy * dy) + t = mod (360. + RADTODEG (atan2 (dy, dx)), 360.) + call fprintf (IE_LOGFD(ie), + " %7.f %7.2f %7.2f %7.2f %7.2f %5.1f") + call pargr (IE_XORIGIN(ie)) + call pargr (IE_YORIGIN(ie)) + call pargr (dx) + call pargr (dy) + call pargr (r) + call pargr (t) + } + call fprintf (IE_LOGFD(ie), "\n") + case 'y': # Set new origin + IE_XORIGIN(ie) = x + IE_YORIGIN(ie) = y + call fprintf (IE_LOGFD(ie), "Origin: %.2f %.2f\n") + call pargr (IE_XORIGIN(ie)) + call pargr (IE_YORIGIN(ie)) + } + } + + # Print in WCS if necessary. + call ie_mwctran (ie, x, y, wx, wy) + if (x == wx && y == wy) + return + call ie_mwctran (ie, IE_XORIGIN(ie), IE_YORIGIN(ie), xo, yo) + + switch (key) { + case 'x': # Print position and pixel value + if (IE_XFORMAT(ie) == '%') + call printf (IE_XFORMAT(ie)) + else + call printf ("%7g") + call pargr (wx) + call printf (" ") + if (IE_YFORMAT(ie) == '%') + call printf (IE_YFORMAT(ie)) + else + call printf ("%7g") + call pargr (wy) + call printf (" %7g") + call pargr (Memr[data]) + + # Print additional fields + if (IE_XORIGIN(ie) != 0. || IE_YORIGIN(ie) != 0.) { + dx = wx - xo + dy = wy - yo + r = sqrt (dx * dx + dy * dy) + t = mod (360. + RADTODEG (atan2 (dy, dx)), 360.) + call printf (" %7g %7g %7g %7g %7g %5.1f") + call pargr (xo) + call pargr (yo) + call pargr (dx) + call pargr (dy) + call pargr (r) + call pargr (t) + } + call printf ("\n") + case 'y': # Set new origin + call printf ("Origin: %7g %7g\n") + call pargr (xo) + call pargr (yo) + } + + # Print to logfile if needed. + if (IE_LOGFD(ie) != NULL) { + switch (key) { + case 'x': + if (IE_XFORMAT(ie) == '%') + call fprintf (IE_LOGFD(ie), IE_XFORMAT(ie)) + else + call fprintf (IE_LOGFD(ie), "%7g") + call pargr (wx) + call fprintf (IE_LOGFD(ie), " ") + if (IE_YFORMAT(ie) == '%') + call fprintf (IE_LOGFD(ie), IE_YFORMAT(ie)) + else + call fprintf (IE_LOGFD(ie), "%7g") + call pargr (wy) + call fprintf (IE_LOGFD(ie), " %7g") + call pargr (Memr[data]) + + if (IE_XORIGIN(ie) != 0. || IE_YORIGIN(ie) != 0.) { + dx = wx - xo + dy = wy - yo + r = sqrt (dx * dx + dy * dy) + t = mod (360. + RADTODEG (atan2 (dy, dx)), 360.) + call fprintf (IE_LOGFD(ie), + " %7g %7g %7g %7g %7g %5.1f") + call pargr (xo) + call pargr (yo) + call pargr (dx) + call pargr (dy) + call pargr (r) + call pargr (t) + } + call fprintf (IE_LOGFD(ie), "\n") + case 'y': # Set new origin + call fprintf (IE_LOGFD(ie), "Origin: %7g %7g\n") + call pargr (xo) + call pargr (yo) + } + } +end diff --git a/pkg/images/tv/imexamine/ieprint.x b/pkg/images/tv/imexamine/ieprint.x new file mode 100644 index 00000000..0a7a7602 --- /dev/null +++ b/pkg/images/tv/imexamine/ieprint.x @@ -0,0 +1,67 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include "imexam.h" + +# IE_PRINT -- Print box of pixel values + +procedure ie_print (ie, x, y) + +pointer ie # IMEXAM structure +real x, y # Center of box + +int i, j, x1, x2, y1, y2, nx +pointer im, data, ie_gimage(), ie_gdata() + +begin + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + x1 = x - 5 + 0.5 + x2 = x + 5 + 0.5 + y1 = y - 5 + 0.5 + y2 = y + 5 + 0.5 + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + return + } + nx = x2 - x1 + 1 + + call printf ("%4w") + do i = x1, x2 { + call printf (" %4d ") + call pargi (i) + } + call printf ("\n") + + do j = y2, y1, -1 { + call printf ("%4d") + call pargi (j) + do i = x1, x2 { + call printf (" %5g") + call pargr (Memr[data+(j-y1)*nx+(i-x1)]) + } + call printf ("\n") + } + + if (IE_LOGFD(ie) != NULL) { + call fprintf (IE_LOGFD(ie), "%4w") + do i = x1, x2 { + call fprintf (IE_LOGFD(ie), " %4d ") + call pargi (i) + } + call fprintf (IE_LOGFD(ie), "\n") + + do j = y2, y1, -1 { + call fprintf (IE_LOGFD(ie), "%4d") + call pargi (j) + do i = x1, x2 { + call fprintf (IE_LOGFD(ie), " %5g") + call pargr (Memr[data+(j-y1)*nx+(i-x1)]) + } + call fprintf (IE_LOGFD(ie), "\n") + } + } +end diff --git a/pkg/images/tv/imexamine/ieqrimexam.x b/pkg/images/tv/imexamine/ieqrimexam.x new file mode 100644 index 00000000..68388874 --- /dev/null +++ b/pkg/images/tv/imexamine/ieqrimexam.x @@ -0,0 +1,489 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include <gset.h> +include <math.h> +include <math/gsurfit.h> +include <math/nlfit.h> +include "imexam.h" + +define FITTYPES "|gaussian|moffat|" +define FITGAUSS 1 +define FITMOFFAT 2 + + +# IE_QRIMEXAM -- Radial profile plot and photometry parameters. +# If no GIO pointer is given then only the photometry parameters are printed. +# First find the center using the marginal distributions. Then subtract +# a fit to the background. Compute the moments within the aperture and +# fit a gaussian of fixed center and zero background. Make the plot +# and print the photometry values. + +procedure ie_qrimexam (gp, mode, ie, x, y) + +pointer gp +pointer ie +int mode +real x, y + +bool center, background, medsky, fitplot, clgpsetb() +real radius, buffer, width, magzero, rplot, beta, clgpsetr() +int fittype, xorder, yorder, clgpseti(), strdic() + +int i, j, ns, no, np, nx, ny, npts, x1, x2, y1, y2 +int plist[3], nplist +real bkg, xcntr, ycntr, mag, e, pa, zcntr, wxcntr, wycntr +real params[3] +real fwhm, dfwhm +pointer sp, fittypes, title, coords, im, data, pp, ws, xs, ys, zs, gs, ptr, nl +double sumo, sums, sumxx, sumyy, sumxy +real r, r1, r2, r3, dx, dy, gseval(), amedr() +pointer clopset(), ie_gimage(), ie_gdata(), locpr() +extern ie_gauss(), ie_dgauss(), ie_moffat(), ie_dmoffat() +errchk nlinit, nlfit + +string glabel "#\ + COL LINE RMAG FLUX SKY N RMOM ELLIP PA PEAK GFWHM\n" +string mlabel "#\ + COL LINE RMAG FLUX SKY N RMOM ELLIP PA PEAK MFWHM\n" + +begin + call smark (sp) + call salloc (fittypes, SZ_FNAME, TY_CHAR) + call salloc (title, IE_SZTITLE, TY_CHAR) + call salloc (coords, IE_SZTITLE, TY_CHAR) + + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + call sfree (sp) + return + } + + # Open parameter set. + if (gp != NULL) { + if (IE_PP(ie) != NULL) + call clcpset (IE_PP(ie)) + } + pp = clopset ("rimexam") + + center = clgpsetb (pp, "center") + background = clgpsetb (pp, "background") + radius = clgpsetr (pp, "radius") + buffer = clgpsetr (pp, "buffer") + width = clgpsetr (pp, "width") + xorder = clgpseti (pp, "xorder") + yorder = clgpseti (pp, "yorder") + medsky = (xorder <= 0 || yorder <= 0) + + magzero = clgpsetr (pp, "magzero") + rplot = clgpsetr (pp, "rplot") + fitplot = clgpsetb (pp, "fitplot") + call clgpseta (pp, "fittype", Memc[fittypes], SZ_FNAME) + fittype = strdic (Memc[fittypes], Memc[fittypes], SZ_FNAME, FITTYPES) + if (fittype == 0) { + call eprintf ("WARNING: Unknown profile fit type `%s'.\n") + call pargstr (Memc[fittypes]) + call sfree (sp) + return + } + beta = clgpsetr (pp, "beta") + + # If the initial center is INDEF then use the previous value. + if (gp != NULL) { + if (!IS_INDEF(x)) + IE_X1(ie) = x + if (!IS_INDEF(y)) + IE_Y1(ie) = y + + xcntr = IE_X1(ie) + ycntr = IE_Y1(ie) + } else { + xcntr = x + ycntr = y + } + + # Center + if (center) + iferr (call ie_center (im, radius, xcntr, ycntr)) { + call erract (EA_WARN) + return + } + + # Crude estimage of FHWM. + dfwhm = radius + + # Get data including a buffer and background annulus. + if (!background) { + buffer = 0. + width = 0. + } + r = max (rplot, radius + buffer + width) + x1 = xcntr - r + x2 = xcntr + r + y1 = ycntr - r + y2 = ycntr + r + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + call sfree (sp) + return + } + + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + npts = nx * ny + + call salloc (xs, npts, TY_REAL) + call salloc (ys, npts, TY_REAL) + call salloc (ws, npts, TY_REAL) + + # Extract the background data if background subtracting. + ns = 0 + if (background && width > 0.) { + call salloc (zs, npts, TY_REAL) + + r1 = radius ** 2 + r2 = (radius + buffer) ** 2 + r3 = (radius + buffer + width) ** 2 + + ptr = data + do j = y1, y2 { + dy = (ycntr - j) ** 2 + do i = x1, x2 { + r = (xcntr - i) ** 2 + dy + if (r <= r1) + ; + else if (r >= r2 && r <= r3) { + Memr[xs+ns] = i + Memr[ys+ns] = j + Memr[zs+ns] = Memr[ptr] + ns = ns + 1 + } + ptr = ptr + 1 + } + } + } + + # Accumulate the various sums for the moments and the gaussian fit. + no = 0 + np = 0 + zcntr = 0. + sumo = 0.; sums = 0.; sumxx = 0.; sumyy = 0.; sumxy = 0. + ptr = data + gs = NULL + + if (ns > 0) { # Background subtraction + + # If background points are defined fit a surface and subtract + # the fitted background from within the object aperture. + + if (medsky) + bkg = amedr (Memr[zs], ns) + else { + repeat { + call gsinit (gs, GS_POLYNOMIAL, xorder, yorder, YES, + real (x1), real (x2), real (y1), real (y2)) + call gsfit (gs, Memr[xs], Memr[ys], Memr[zs], Memr[ws], ns, + WTS_UNIFORM, i) + if (i == OK) + break + xorder = max (1, xorder - 1) + yorder = max (1, yorder - 1) + call gsfree (gs) + } + bkg = gseval (gs, real(x1), real(y1)) + } + + do j = y1, y2 { + dy = j - ycntr + do i = x1, x2 { + dx = i - xcntr + r = sqrt (dx ** 2 + dy ** 2) + r3 = max (0., min (5., 2 * r / dfwhm - 1.)) + + if (medsky) + r2 = bkg + else { + r2 = gseval (gs, real(i), real(j)) + bkg = min (bkg, r2) + } + r1 = Memr[ptr] - r2 + + if (r <= radius) { + sumo = sumo + r1 + sums = sums + r2 + sumxx = sumxx + dx * dx * r1 + sumyy = sumyy + dy * dy * r1 + sumxy = sumxy + dx * dy * r1 + zcntr = max (r1, zcntr) + if (r <= rplot) { + Memr[xs+no] = r + Memr[ys+no] = r1 + Memr[ws+no] = exp (-r3**2) / max (.1, r**2) + no = no + 1 + } else { + np = np + 1 + Memr[xs+npts-np] = r + Memr[ys+npts-np] = r1 + Memr[ws+npts-np] = exp (-r3**2) / max (.1, r**2) + } + } else if (r <= rplot) { + np = np + 1 + Memr[xs+npts-np] = r + Memr[ys+npts-np] = r1 + } + ptr = ptr + 1 + } + } + + if (gs != NULL) + call gsfree (gs) + + } else { # No background subtraction + bkg = 0. + do j = y1, y2 { + dy = j - ycntr + do i = x1, x2 { + dx = i - xcntr + r = sqrt (dx ** 2 + dy ** 2) + r3 = max (0., min (5., 2 * r / dfwhm - 1.)) + r1 = Memr[ptr] + + if (r <= radius) { + sumo = sumo + r1 + sumxx = sumxx + dx * dx * r1 + sumyy = sumyy + dy * dy * r1 + sumxy = sumxy + dx * dy * r1 + zcntr = max (r1, zcntr) + if (r <= rplot) { + Memr[xs+no] = r + Memr[ys+no] = r1 + Memr[ws+no] = exp (-r3**2) / max (.1, r**2) + no = no + 1 + } else { + np = np + 1 + Memr[xs+npts-np] = r + Memr[ys+npts-np] = r1 + Memr[ws+npts-np] = exp (-r3**2) / max (.1, r**2) + } + } else if (r <= rplot) { + np = np + 1 + Memr[xs+npts-np] = r + Memr[ys+npts-np] = r1 + } + ptr = ptr + 1 + } + } + } + if (np > 0) { + call amovr (Memr[xs+npts-np], Memr[xs+no], np) + call amovr (Memr[ys+npts-np], Memr[ys+no], np) + call amovr (Memr[ws+npts-np], Memr[ws+no], np) + } + if (rplot <= radius) { + no = no + np + np = no - np + } else + np = no + np + + + # Compute the photometry and gaussian fit parameters. + + switch (fittype) { + case FITGAUSS: + plist[1] = 1 + plist[2] = 2 + nplist = 2 + params[2] = dfwhm**2 / (8 * log(2.)) + params[1] = zcntr + call nlinitr (nl, locpr (ie_gauss), locpr (ie_dgauss), + params, params, 2, plist, nplist, .001, 100) + call nlfitr (nl, Memr[xs], Memr[ys], Memr[ws], no, 1, WTS_USER, i) + if (i == SINGULAR || i == NO_DEG_FREEDOM) { + call eprintf ("WARNING: Gaussian fit did not converge\n") + call tsleep (5) + zcntr = INDEF + fwhm = INDEF + } else { + call nlpgetr (nl, params, i) + if (params[2] < 0.) { + zcntr = INDEF + fwhm = INDEF + } else { + zcntr = params[1] + fwhm = sqrt (8 * log (2.) * params[2]) + } + } + case FITMOFFAT: + plist[1] = 1 + plist[2] = 2 + if (IS_INDEF(beta)) { + params[3] = -3.0 + plist[3] = 3 + nplist = 3 + } else { + params[3] = -beta + nplist = 2 + } + params[2] = dfwhm / 2. / sqrt (2.**(-1./params[3]) - 1.) + params[1] = zcntr + call nlinitr (nl, locpr (ie_moffat), locpr (ie_dmoffat), + params, params, 3, plist, nplist, .001, 100) + call nlfitr (nl, Memr[xs], Memr[ys], Memr[ws], no, 1, WTS_USER, i) + if (i == SINGULAR || i == NO_DEG_FREEDOM) { + call eprintf ("WARNING: Moffat fit did not converge\n") + call tsleep (5) + zcntr = INDEF + fwhm = INDEF + beta = INDEF + } else { + call nlpgetr (nl, params, i) + if (params[2] < 0.) { + zcntr = INDEF + fwhm = INDEF + beta = INDEF + } else { + zcntr = params[1] + beta = -params[3] + fwhm = abs (params[2])*2.*sqrt (2.**(-1./params[3]) - 1.) + } + } + } + + mag = INDEF + r = INDEF + e = INDEF + pa = INDEF + if (sumo > 0.) { + mag = magzero - 2.5 * log10 (sumo) + r2 = sumxx + sumyy + if (r2 > 0.) { + switch (fittype) { + case FITGAUSS: + r = 2 * sqrt (log (2.) * r2 / sumo) + case FITMOFFAT: + if (beta > 2.) + r = 2 * sqrt ((beta-2.)*(2.**(1./beta)-1) * r2 / sumo) + } + r1 =(sumxx-sumyy)**2+(2*sumxy)**2 + if (r1 > 0.) + e = sqrt (r1) / r2 + else + e = 0. + } + if (e < 0.01) + e = 0. + else + pa = RADTODEG (0.5 * atan2 (2*sumxy, sumxx-sumyy)) + } + + call ie_mwctran (ie, xcntr, ycntr, wxcntr, wycntr) + if (xcntr == wxcntr && ycntr == wycntr) + call strcpy ("%.2f %.2f", Memc[title], IE_SZTITLE) + else { + call sprintf (Memc[title], IE_SZTITLE, "%s %s") + if (IE_XFORMAT(ie) == '%') + call pargstr (IE_XFORMAT(ie)) + else + call pargstr ("%g") + if (IE_YFORMAT(ie) == '%') + call pargstr (IE_YFORMAT(ie)) + else + call pargstr ("%g") + } + call sprintf (Memc[coords], IE_SZTITLE, Memc[title]) + call pargr (wxcntr) + call pargr (wycntr) + + # Plot the radial profile and overplot the fit. + if (gp != NULL) { + call sprintf (Memc[title], IE_SZTITLE, + "%s: Radial profile at %s\n%s") + call pargstr (IE_IMNAME(ie)) + call pargstr (Memc[coords]) + call pargstr (IM_TITLE(im)) + + call ie_graph (gp, mode, pp, Memc[title], Memr[xs], Memr[ys], + np, "", "") + + if (fitplot && !IS_INDEF (fwhm)) { + np = 51 + dx = rplot / (np - 1) + do i = 0, np - 1 + Memr[xs+i] = i * dx + call nlvectorr (nl, Memr[xs], Memr[ys], np, 1) + call gseti (gp, G_PLTYPE, 2) + call gpline (gp, Memr[xs], Memr[ys], np) + call gseti (gp, G_PLTYPE, 1) + } + } + + if (IE_LASTKEY(ie) != ',') { + switch (fittype) { + case FITGAUSS: + call printf (glabel) + case FITMOFFAT: + call printf (mlabel) + } + } + + # Print the photometry values. + call printf ( + "%7.2f %7.2f %7.2f %8.1f %8.2f %3d %5.2f %5.3f %5.1f %8.2f %5.2f\n") + call pargr (xcntr) + call pargr (ycntr) + call pargr (mag) + call pargd (sumo) + call pargd (sums / no) + call pargi (no) + call pargr (r) + call pargr (e) + call pargr (pa) + call pargr (zcntr) + call pargr (fwhm) + if (gp == NULL) { + if (xcntr != wxcntr || ycntr != wycntr) { + call printf ("%s: %s\n") + call pargstr (IE_WCSNAME(ie)) + call pargstr (Memc[coords]) + } + } + + if (IE_LOGFD(ie) != NULL) { + if (IE_LASTKEY(ie) != ',') { + switch (fittype) { + case FITGAUSS: + call fprintf (IE_LOGFD(ie), glabel) + case FITMOFFAT: + call fprintf (IE_LOGFD(ie), mlabel) + } + } + + call fprintf (IE_LOGFD(ie), + "%7.2f %7.2f %7.2f %8.1f %8.2f %3d %5.2f %5.3f %5.1f %8.2f %5.2f\n") + call pargr (xcntr) + call pargr (ycntr) + call pargr (mag) + call pargd (sumo) + call pargd (sums / no) + call pargi (no) + call pargr (r) + call pargr (e) + call pargr (pa) + call pargr (zcntr) + call pargr (fwhm) + if (xcntr != wxcntr || ycntr != wycntr) { + call fprintf (IE_LOGFD(ie), "%s: %s\n") + call pargstr (IE_WCSNAME(ie)) + call pargstr (Memc[coords]) + } + } + + if (gp == NULL) + call clcpset (pp) + else + IE_PP(ie) = pp + + call nlfreer (nl) + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/ierimexam.x b/pkg/images/tv/imexamine/ierimexam.x new file mode 100644 index 00000000..f76ff507 --- /dev/null +++ b/pkg/images/tv/imexamine/ierimexam.x @@ -0,0 +1,752 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include <gset.h> +include <math.h> +include <math/gsurfit.h> +include <math/nlfit.h> +include "imexam.h" + +define FITTYPES "|gaussian|moffat|" +define FITGAUSS 1 +define FITMOFFAT 2 + + +# IE_RIMEXAM -- Radial profile plot and photometry parameters. +# If no GIO pointer is given then only the photometry parameters are printed. +# First find the center using the marginal distributions. Then subtract +# a fit to the background. Compute the moments within the aperture and +# fit a gaussian of fixed center and zero background. Make the plot +# and print the photometry values. + +procedure ie_rimexam (gp, mode, ie, x, y) + +pointer gp +pointer ie +int mode +real x, y + +bool center, background, medsky, fitplot, clgpsetb() +real radius, buffer, width, magzero, rplot, beta, clgpsetr() +int nit, fittype, xorder, yorder, clgpseti(), strdic() + +int i, j, ns, no, np, nx, ny, npts, x1, x2, y1, y2 +int coordlen, plist[3], nplist, strlen() +real bkg, xcntr, ycntr, mag, e, pa, zcntr, wxcntr, wycntr +real params[3] +real fwhm, dbkg, dfwhm, gfwhm, efwhm +pointer sp, fittypes, title, coords, im, data, pp, ws, xs, ys, zs, gs, ptr, nl +double sumo, sums, sumxx, sumyy, sumxy +real r, r1, r2, r3, dx, dy, gseval(), amedr() +pointer clopset(), ie_gimage(), ie_gdata(), locpr() +extern ie_gauss(), ie_dgauss(), ie_moffat(), ie_dmoffat() +errchk stf_measure, nlinit, nlfit + +begin + call smark (sp) + call salloc (fittypes, SZ_FNAME, TY_CHAR) + call salloc (title, IE_SZTITLE, TY_CHAR) + call salloc (coords, IE_SZTITLE, TY_CHAR) + + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + call sfree (sp) + return + } + + # Open parameter set. + if (gp != NULL) { + if (IE_PP(ie) != NULL) + call clcpset (IE_PP(ie)) + } + pp = clopset ("rimexam") + + center = clgpsetb (pp, "center") + background = clgpsetb (pp, "background") + radius = clgpsetr (pp, "radius") + buffer = clgpsetr (pp, "buffer") + width = clgpsetr (pp, "width") + xorder = clgpseti (pp, "xorder") + yorder = clgpseti (pp, "yorder") + medsky = (xorder <= 0 || yorder <= 0) + nit = clgpseti (pp, "iterations") + + magzero = clgpsetr (pp, "magzero") + rplot = clgpsetr (pp, "rplot") + fitplot = clgpsetb (pp, "fitplot") + call clgpseta (pp, "fittype", Memc[fittypes], SZ_FNAME) + fittype = strdic (Memc[fittypes], Memc[fittypes], SZ_FNAME, FITTYPES) + if (fittype == 0) { + call eprintf ("WARNING: Unknown profile fit type `%s'.\n") + call pargstr (Memc[fittypes]) + call sfree (sp) + return + } + beta = clgpsetr (pp, "beta") + + # If the initial center is INDEF then use the previous value. + if (gp != NULL) { + if (!IS_INDEF(x)) + IE_X1(ie) = x + if (!IS_INDEF(y)) + IE_Y1(ie) = y + + xcntr = IE_X1(ie) + ycntr = IE_Y1(ie) + } else { + xcntr = x + ycntr = y + } + + # Center + if (center) + iferr (call ie_center (im, radius, xcntr, ycntr)) { + call erract (EA_WARN) + call sfree (sp) + return + } + + # Do the enclosed flux and direct FWHM measurments using the + # PSFMEASURE routines. + + call stf_measure (im, xcntr, ycntr, beta, 0.5, radius, nit, buffer, + width, INDEF, NULL, NULL, dbkg, r, dfwhm, gfwhm, efwhm) + if (fittype == FITGAUSS) + efwhm = gfwhm + + # Get data including a buffer and background annulus. + if (!background) { + buffer = 0. + width = 0. + } + r = max (rplot, radius + buffer + width) + x1 = xcntr - r + x2 = xcntr + r + y1 = ycntr - r + y2 = ycntr + r + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + call sfree (sp) + return + } + + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + npts = nx * ny + + call salloc (xs, npts, TY_REAL) + call salloc (ys, npts, TY_REAL) + call salloc (ws, npts, TY_REAL) + + # Extract the background data if background subtracting. + ns = 0 + if (background && width > 0.) { + call salloc (zs, npts, TY_REAL) + + r1 = radius ** 2 + r2 = (radius + buffer) ** 2 + r3 = (radius + buffer + width) ** 2 + + ptr = data + do j = y1, y2 { + dy = (ycntr - j) ** 2 + do i = x1, x2 { + r = (xcntr - i) ** 2 + dy + if (r <= r1) + ; + else if (r >= r2 && r <= r3) { + Memr[xs+ns] = i + Memr[ys+ns] = j + Memr[zs+ns] = Memr[ptr] + ns = ns + 1 + } + ptr = ptr + 1 + } + } + } + + # Accumulate the various sums for the moments and the gaussian fit. + no = 0 + np = 0 + zcntr = 0. + sumo = 0.; sums = 0.; sumxx = 0.; sumyy = 0.; sumxy = 0. + ptr = data + gs = NULL + + if (ns > 0) { # Background subtraction + + # If background points are defined fit a surface and subtract + # the fitted background from within the object aperture. + + if (medsky) + bkg = amedr (Memr[zs], ns) + else { + repeat { + call gsinit (gs, GS_POLYNOMIAL, xorder, yorder, YES, + real (x1), real (x2), real (y1), real (y2)) + call gsfit (gs, Memr[xs], Memr[ys], Memr[zs], Memr[ws], ns, + WTS_UNIFORM, i) + if (i == OK) + break + xorder = max (1, xorder - 1) + yorder = max (1, yorder - 1) + call gsfree (gs) + } + bkg = gseval (gs, real(x1), real(y1)) + } + + do j = y1, y2 { + dy = j - ycntr + do i = x1, x2 { + dx = i - xcntr + r = sqrt (dx ** 2 + dy ** 2) + r3 = max (0., min (5., 2 * r / dfwhm - 1.)) + + if (medsky) + r2 = bkg + else { + r2 = gseval (gs, real(i), real(j)) + bkg = min (bkg, r2) + } + r1 = Memr[ptr] - r2 + + if (r <= radius) { + sumo = sumo + r1 + sums = sums + r2 + sumxx = sumxx + dx * dx * r1 + sumyy = sumyy + dy * dy * r1 + sumxy = sumxy + dx * dy * r1 + zcntr = max (r1, zcntr) + if (r <= rplot) { + Memr[xs+no] = r + Memr[ys+no] = r1 + Memr[ws+no] = exp (-r3**2) / max (.1, r**2) + no = no + 1 + } else { + np = np + 1 + Memr[xs+npts-np] = r + Memr[ys+npts-np] = r1 + Memr[ws+npts-np] = exp (-r3**2) / max (.1, r**2) + } + } else if (r <= rplot) { + np = np + 1 + Memr[xs+npts-np] = r + Memr[ys+npts-np] = r1 + } + ptr = ptr + 1 + } + } + + if (gs != NULL) + call gsfree (gs) + + } else { # No background subtraction + bkg = 0. + do j = y1, y2 { + dy = j - ycntr + do i = x1, x2 { + dx = i - xcntr + r = sqrt (dx ** 2 + dy ** 2) + r3 = max (0., min (5., 2 * r / dfwhm - 1.)) + r1 = Memr[ptr] + + if (r <= radius) { + sumo = sumo + r1 + sumxx = sumxx + dx * dx * r1 + sumyy = sumyy + dy * dy * r1 + sumxy = sumxy + dx * dy * r1 + zcntr = max (r1, zcntr) + if (r <= rplot) { + Memr[xs+no] = r + Memr[ys+no] = r1 + Memr[ws+no] = exp (-r3**2) / max (.1, r**2) + no = no + 1 + } else { + np = np + 1 + Memr[xs+npts-np] = r + Memr[ys+npts-np] = r1 + Memr[ws+npts-np] = exp (-r3**2) / max (.1, r**2) + } + } else if (r <= rplot) { + np = np + 1 + Memr[xs+npts-np] = r + Memr[ys+npts-np] = r1 + } + ptr = ptr + 1 + } + } + } + if (np > 0) { + call amovr (Memr[xs+npts-np], Memr[xs+no], np) + call amovr (Memr[ys+npts-np], Memr[ys+no], np) + call amovr (Memr[ws+npts-np], Memr[ws+no], np) + } + if (rplot <= radius) { + no = no + np + np = no - np + } else + np = no + np + + + # Compute the photometry and profile fit parameters. + + switch (fittype) { + case FITGAUSS: + plist[1] = 1 + plist[2] = 2 + nplist = 2 + params[2] = dfwhm**2 / (8 * log(2.)) + params[1] = zcntr + call nlinitr (nl, locpr (ie_gauss), locpr (ie_dgauss), + params, params, 2, plist, nplist, .001, 100) + call nlfitr (nl, Memr[xs], Memr[ys], Memr[ws], no, 1, WTS_USER, i) + if (i == SINGULAR || i == NO_DEG_FREEDOM) { + call eprintf ("WARNING: Gaussian fit did not converge\n") + call tsleep (5) + zcntr = INDEF + fwhm = INDEF + } else { + call nlpgetr (nl, params, i) + if (params[2] < 0.) { + zcntr = INDEF + fwhm = INDEF + } else { + zcntr = params[1] + fwhm = sqrt (8 * log (2.) * params[2]) + } + } + case FITMOFFAT: + plist[1] = 1 + plist[2] = 2 + if (IS_INDEF(beta)) { + params[3] = -3.0 + plist[3] = 3 + nplist = 3 + } else { + params[3] = -beta + nplist = 2 + } + params[2] = dfwhm / 2. / sqrt (2.**(-1./params[3]) - 1.) + params[1] = zcntr + call nlinitr (nl, locpr (ie_moffat), locpr (ie_dmoffat), + params, params, 3, plist, nplist, .001, 100) + call nlfitr (nl, Memr[xs], Memr[ys], Memr[ws], no, 1, WTS_USER, i) + if (i == SINGULAR || i == NO_DEG_FREEDOM) { + call eprintf ("WARNING: Moffat fit did not converge\n") + call tsleep (5) + zcntr = INDEF + fwhm = INDEF + beta = INDEF + } else { + call nlpgetr (nl, params, i) + if (params[2] < 0.) { + zcntr = INDEF + fwhm = INDEF + beta = INDEF + } else { + zcntr = params[1] + beta = -params[3] + fwhm = abs (params[2])*2.*sqrt (2.**(-1./params[3]) - 1.) + } + } + } + + mag = INDEF + r = INDEF + e = INDEF + pa = INDEF + if (sumo > 0.) { + mag = magzero - 2.5 * log10 (sumo) + r2 = sumxx + sumyy + if (r2 > 0.) { + switch (fittype) { + case FITGAUSS: + r = 2 * sqrt (log (2.) * r2 / sumo) + case FITMOFFAT: + if (beta > 2.) + r = 2 * sqrt ((beta-2.)*(2.**(1./beta)-1) * r2 / sumo) + } + r1 =(sumxx-sumyy)**2+(2*sumxy)**2 + if (r1 > 0.) + e = sqrt (r1) / r2 + else + e = 0. + } + if (e < 0.01) + e = 0. + else + pa = RADTODEG (0.5 * atan2 (2*sumxy, sumxx-sumyy)) + } + + call ie_mwctran (ie, xcntr, ycntr, wxcntr, wycntr) + if (xcntr == wxcntr && ycntr == wycntr) + call strcpy ("%.2f %.2f", Memc[title], IE_SZTITLE) + else { + call sprintf (Memc[title], IE_SZTITLE, "%s %s") + if (IE_XFORMAT(ie) == '%') + call pargstr (IE_XFORMAT(ie)) + else + call pargstr ("%g") + if (IE_YFORMAT(ie) == '%') + call pargstr (IE_YFORMAT(ie)) + else + call pargstr ("%g") + } + call sprintf (Memc[coords], IE_SZTITLE, Memc[title]) + call pargr (wxcntr) + call pargr (wycntr) + + # Plot the radial profile and overplot the gaussian fit. + if (gp != NULL) { + call sprintf (Memc[title], IE_SZTITLE, + "%s: Radial profile at %s\n%s") + call pargstr (IE_IMNAME(ie)) + call pargstr (Memc[coords]) + call pargstr (IM_TITLE(im)) + + call ie_graph (gp, mode, pp, Memc[title], Memr[xs], Memr[ys], + np, "", "") + + if (fitplot && !IS_INDEF (fwhm)) { + np = 51 + dx = rplot / (np - 1) + do i = 0, np - 1 + Memr[xs+i] = i * dx + call nlvectorr (nl, Memr[xs], Memr[ys], np, 1) + call gseti (gp, G_PLTYPE, 2) + call gpline (gp, Memr[xs], Memr[ys], np) + call gseti (gp, G_PLTYPE, 1) + } + call gseti (gp, G_PLTYPE, 2) + + call printf ("%6.2f %6.2f %7.4g %7.4g %7.4g %4.2f %4d") + call pargr (radius) + call pargr (mag) + call pargd (sumo) + call pargd (sums / no) + call pargr (zcntr) + call pargr (e) + call pargr (pa) + switch (fittype) { + case FITGAUSS: + call printf (" %4w %8.2f %8.2f %6.2f\n") + call pargr (efwhm) + call pargr (fwhm) + call pargr (dfwhm) + case FITMOFFAT: + call printf (" %4.2f %8.2f %8.2f %6.2f\n") + call pargr (beta) + call pargr (efwhm) + call pargr (fwhm) + call pargr (dfwhm) + } + + } else { + if (IE_LASTKEY(ie) != 'a') { + coordlen = max (11, strlen (Memc[coords])) + call printf ("# %5s %7s %-*s\n# %5s %6s %7s %7s %7s %4s %4s") + call pargstr ("COL") + call pargstr ("LINE") + call pargi (coordlen) + call pargstr ("COORDINATES") + call pargstr ("R") + call pargstr ("MAG") + call pargstr ("FLUX") + call pargstr ("SKY") + call pargstr ("PEAK") + call pargstr ("E") + call pargstr ("PA") + switch (fittype) { + case FITGAUSS: + call printf (" %4w %8s %8s %6s\n") + call pargstr ("ENCLOSED") + call pargstr ("GAUSSIAN") + call pargstr ("DIRECT") + case FITMOFFAT: + call printf (" %4s %8s %8s %6s\n") + call pargstr ("BETA") + call pargstr ("ENCLOSED") + call pargstr ("MOFFAT") + call pargstr ("DIRECT") + } + } + + call printf ( + "%7.2f %7.2f %-*s\n %6.2f %6.2f %7.4g %7.4g %7.4g %4.2f %4d") + call pargr (xcntr) + call pargr (ycntr) + call pargi (coordlen) + call pargstr (Memc[coords]) + call pargr (radius) + call pargr (mag) + call pargd (sumo) + call pargd (sums / no) + call pargr (zcntr) + call pargr (e) + call pargr (pa) + switch (fittype) { + case FITGAUSS: + call printf (" %4w %8.2f %8.2f %6.2f\n") + call pargr (efwhm) + call pargr (fwhm) + call pargr (dfwhm) + case FITMOFFAT: + call printf (" %4.2f %8.2f %8.2f %6.2f\n") + call pargr (beta) + call pargr (efwhm) + call pargr (fwhm) + call pargr (dfwhm) + } + } + + if (IE_LOGFD(ie) != NULL) { + if (IE_LASTKEY(ie) != 'a') { + coordlen = max (11, strlen (Memc[coords])) + call fprintf (IE_LOGFD(ie), + "# %5s %7s %-*s %6s %6s %7s %7s %7s %4s %4s") + call pargstr ("COL") + call pargstr ("LINE") + call pargi (coordlen) + call pargstr ("COORDINATES") + call pargstr ("R") + call pargstr ("MAG") + call pargstr ("FLUX") + call pargstr ("SKY") + call pargstr ("PEAK") + call pargstr ("E") + call pargstr ("PA") + switch (fittype) { + case FITGAUSS: + call fprintf (IE_LOGFD(ie), " %4w %8s %8s %6s\n") + call pargstr ("ENCLOSED") + call pargstr ("GAUSSIAN") + call pargstr ("DIRECT") + case FITMOFFAT: + call fprintf (IE_LOGFD(ie), " %4s %8s %8s %6s\n") + call pargstr ("BETA") + call pargstr ("ENCLOSED") + call pargstr ("MOFFAT") + call pargstr ("DIRECT") + } + } + + call fprintf (IE_LOGFD(ie), + "%7.2f %7.2f %-*s %6.2f %6.2f %7.4g %7.4g %7.4g %4.2f %4d") + call pargr (xcntr) + call pargr (ycntr) + call pargi (coordlen) + call pargstr (Memc[coords]) + call pargr (radius) + call pargr (mag) + call pargd (sumo) + call pargd (sums / no) + call pargr (zcntr) + call pargr (e) + call pargr (pa) + switch (fittype) { + case FITGAUSS: + call fprintf (IE_LOGFD(ie), " %4w %8.2f %8.2f %6.2f\n") + call pargr (efwhm) + call pargr (fwhm) + call pargr (dfwhm) + case FITMOFFAT: + call fprintf (IE_LOGFD(ie), " %4.2f %8.2f %8.2f %6.2f\n") + call pargr (beta) + call pargr (efwhm) + call pargr (fwhm) + call pargr (dfwhm) + } + } + + if (gp == NULL) + call clcpset (pp) + else + IE_PP(ie) = pp + + call nlfreer (nl) + call sfree (sp) +end + + +# IE_CENTER -- Find the center of gravity from the marginal distributions. + +procedure ie_center (im, radius, xcntr, ycntr) + +pointer im +real radius +real xcntr, ycntr + +int i, j, k, x1, x2, y1, y2, nx, ny, npts +real xlast, ylast +real mean, sum, sum1, sum2, sum3, asumr() +pointer data, ptr, ie_gdata() +errchk ie_gdata + +begin + # Find the center of a star image given approximate coords. Uses + # Mountain Photometry Code Algorithm as outlined in Stellar Magnitudes + # from Digital Images. + + do k = 1, 3 { + # Extract region around center + xlast = xcntr + ylast = ycntr + x1 = xcntr - radius + 0.5 + x2 = xcntr + radius + 0.5 + y1 = ycntr - radius + 0.5 + y2 = ycntr + radius + 0.5 + data = ie_gdata (im, x1, x2, y1, y2) + + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + npts = nx * ny + + # Find center of gravity for marginal distributions above mean. + sum = asumr (Memr[data], npts) + mean = sum / nx + sum1 = 0. + sum2 = 0. + + do i = x1, x2 { + ptr = data + i - x1 + sum3 = 0. + do j = y1, y2 { + sum3 = sum3 + Memr[ptr] + ptr = ptr + nx + } + sum3 = sum3 - mean + if (sum3 > 0.) { + sum1 = sum1 + i * sum3 + sum2 = sum2 + sum3 + } + } + xcntr = sum1 / sum2 + + ptr = data + mean = sum / ny + sum1 = 0. + sum2 = 0. + do j = y1, y2 { + sum3 = 0. + do i = x1, x2 { + sum3 = sum3 + Memr[ptr] + ptr = ptr + 1 + } + sum3 = sum3 - mean + if (sum3 > 0.) { + sum1 = sum1 + j * sum3 + sum2 = sum2 + sum3 + } + } + ycntr = sum1 / sum2 + + if (int(xcntr) == int(xlast) && int(ycntr) == int(ylast)) + break + } +end + + +# IE_GAUSS -- Gaussian function used in NLFIT. The parameters are the +# amplitude and sigma squared and the input variable is the radius. + +procedure ie_gauss (x, nvars, p, np, z) + +real x[nvars] #I Input variables +int nvars #I Number of variables +real p[np] #I Parameter vector +int np #I Number of parameters +real z #O Function return + +real r2 + +begin + r2 = x[1]**2 / (2 * p[2]) + if (abs (r2) > 20.) + z = 0. + else + z = p[1] * exp (-r2) +end + + +# IE_DGAUSS -- Gaussian function and derivatives used in NLFIT. The parameters +# are the amplitude and sigma squared and the input variable is the radius. + +procedure ie_dgauss (x, nvars, p, dp, np, z, der) + +real x[nvars] #I Input variables +int nvars #I Number of variables +real p[np] #I Parameter vector +real dp[np] #I Dummy array of parameters increments +int np #I Number of parameters +real z #O Function return +real der[np] #O Derivatives + +real r2 + +begin + r2 = x[1]**2 / (2 * p[2]) + if (abs (r2) > 20.) { + z = 0. + der[1] = 0. + der[2] = 0. + } else { + der[1] = exp (-r2) + z = p[1] * der[1] + der[2] = z * r2 / p[2] + } +end + + +# IE_MOFFAT -- Moffat function used in NLFIT. The parameters are the +# amplitude, alpha squared, and beta and the input variable is the radius. + +procedure ie_moffat (x, nvars, p, np, z) + +real x[nvars] #I Input variables +int nvars #I Number of variables +real p[np] #I Parameter vector +int np #I Number of parameters +real z #O Function return + +real y + +begin + y = 1 + (x[1] / p[2]) ** 2 + if (abs (y) > 20.) + z = 0. + else + z = p[1] * y ** p[3] +end + + +# IE_DMOFFAT -- Moffat function and derivatives used in NLFIT. The parameters +# are the amplitude, alpha squared, and beta and the input variable is the +# radius. + +procedure ie_dmoffat (x, nvars, p, dp, np, z, der) + +real x[nvars] #I Input variables +int nvars #I Number of variables +real p[np] #I Parameter vector +real dp[np] #I Dummy array of parameters increments +int np #I Number of parameters +real z #O Function return +real der[np] #O Derivatives + +real y + +begin + y = 1 + (x[1] / p[2]) ** 2 + if (abs (y) > 20.) { + z = 0. + der[1] = 0. + der[2] = 0. + der[3] = 0. + } else { + der[1] = y ** p[3] + z = p[1] * der[1] + der[2] = -2 * z / y * p[3] / p[2] * (x[1] / p[2]) ** 2 + der[3] = z * log (y) + } +end diff --git a/pkg/images/tv/imexamine/iesimexam.x b/pkg/images/tv/imexamine/iesimexam.x new file mode 100644 index 00000000..292364ee --- /dev/null +++ b/pkg/images/tv/imexamine/iesimexam.x @@ -0,0 +1,492 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include <gset.h> +include <mach.h> +include "imexam.h" + +define CSIZE 24 + + +# IE_SIMEXAM -- Draw a perspective view of a surface. The altitude +# and azimuth of the viewing angle are variable. + +procedure ie_simexam (gp, mode, ie, x, y) + +pointer gp # GIO pointer +int mode # Mode +pointer ie # IMEXAM pointer +real x, y # Center + +real angh, angv # Orientation of surface (degrees) +real floor, ceiling # Range limits + +int wkid +int x1, x2, y1, y2, nx, ny, npts +pointer pp, sp, title, str, sdata, work, im, data, ie_gimage(), ie_gdata() + +bool clgpsetb() +int clgpseti() +real clgpsetr() +pointer clopset() + +int first +real vpx1, vpx2, vpy1, vpy2 +common /frstfg/ first +common /noaovp/ vpx1, vpx2, vpy1, vpy2 + +begin + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + pp = IE_PP(ie) + if (pp != NULL) + call clcpset (pp) + pp = clopset ("simexam") + IE_PP(ie) = pp + + nx = clgpseti (pp, "ncolumns") + ny = clgpseti (pp, "nlines") + angh = clgpsetr (pp, "angh") + angv = clgpsetr (pp, "angv") + floor = clgpsetr (pp, "floor") + ceiling = clgpsetr (pp, "ceiling") + + if (!IS_INDEF(x)) + IE_X1(ie) = x + if (!IS_INDEF(y)) + IE_Y1(ie) = y + + x1 = IE_X1(ie) - (nx - 1) / 2 + 0.5 + x2 = IE_X1(ie) + nx / 2 + 0.5 + y1 = IE_Y1(ie) - (ny - 1) / 2 + 0.5 + y2 = IE_Y1(ie) + ny / 2 + 0.5 + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + return + } + + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + npts = nx * ny + + call smark (sp) + + # Take floor and ceiling if enabled (nonzero). + if (IS_INDEF (floor) && IS_INDEF (ceiling)) + sdata = data + else { + call salloc (sdata, npts, TY_REAL) + call amovr (Memr[data], Memr[sdata], npts) + if (!IS_INDEF (floor) && !IS_INDEF (ceiling)) { + floor = min (floor, ceiling) + ceiling = max (floor, ceiling) + } + } + iferr (call ie_surf_limits (Memr[sdata], npts, floor, ceiling)) { + call sfree (sp) + call erract (EA_WARN) + return + } + + if (mode != APPEND) { + call gclear (gp) + + # Set the viewport. + call gsview (gp, 0.1, 0.9, 0.1, 0.9) + + call salloc (title, IE_SZTITLE, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + if (clgpsetb (pp, "banner")) { + call sysid (Memc[str], SZ_LINE) + call sprintf (Memc[title], IE_SZTITLE, + "%s\n%s: Surface plot of [%d:%d,%d:%d]\n%s") + call pargstr (Memc[str]) + call pargstr (IE_IMNAME(ie)) + call pargi (x1) + call pargi (x2) + call pargi (y1) + call pargi (y2) + call pargstr (IM_TITLE(im)) + } else + Memc[title] = EOS + + call clgpset (pp, "title", Memc[str], SZ_LINE) + if (Memc[str] != EOS) { + call strcat ("\n", Memc[title], IE_SZTITLE) + call strcat (Memc[str], Memc[title], IE_SZTITLE) + } + + call gseti (gp, G_DRAWAXES, NO) + call glabax (gp, Memc[title], "", "") + } + + # Open graphics device and make plot. + call gopks (STDERR) + wkid = 1 + call gopwk (wkid, 6, gp) + call gacwk (wkid) + + first = 1 + call srfabd() + call ggview (gp, vpx1, vpx2, vpy1, vpy2) + call set (vpx1, vpx2, vpy1, vpy2, 1.0, 1024., 1.0, 1024., 1) + call salloc (work, 2 * (2*nx*ny+nx+ny), TY_REAL) + call ezsrfc (Memr[sdata], nx, ny, angh, angv, Memr[work]) + + if (mode != APPEND) { + if (clgpsetb (pp, "axes")) { + call gswind (gp, real (x1), real (x2), real (y1), real (y2)) + call gseti (gp, G_CLIP, NO) + call ie_perimeter (gp, Memr[sdata], nx, ny, angh, angv) + } + } + + call gdawk (wkid) + call gclks () + call sfree (sp) +end + + +# IE_PERIMETER -- draw and label axes around the surface plot. + +procedure ie_perimeter (gp, z, ncols, nlines, angh, angv) + +pointer gp # Graphics pointer +int ncols # Number of image columns +int nlines # Number of image lines +real z[ncols, nlines] # Array of intensity values +real angh # Angle of horizontal inclination +real angv # Angle of vertical inclination + +pointer sp, x_val, y_val, kvec +char tlabel[10] +real xmin, ymin, delta, fact1, flo, hi, xcen, ycen +real x1_perim, x2_perim, y1_perim, y2_perim, z1, z2 +real wc1, wc2, wl1, wl2, del +int i, j, junk +int itoc() +data fact1 /2.0/ +real vpx1, vpx2, vpy1, vpy2 +common /noaovp/ vpx1, vpx2, vpy1, vpy2 + +begin + call smark (sp) + call salloc (x_val, ncols + 2, TY_REAL) + call salloc (y_val, nlines + 2, TY_REAL) + call salloc (kvec, max (ncols, nlines) + 2, TY_REAL) + + # Get window coordinates set up in calling procedure. + call ggwind (gp, wc1, wc2, wl1, wl2) + + # Set up window, viewport for output. The coordinates returned + # from trn32s are in the range [1-1024]. + call set (vpx1, vpx2, vpy1, vpy2, 1.0, 1024., 1.0, 1024., 1) + + # Find range of z for determining perspective + flo = MAX_REAL + hi = -flo + do j = 1, nlines { + call alimr (z[1,j], ncols, z1, z2) + flo = min (flo, z1) + hi = max (hi, z2) + } + + # Set up linear endpoints and spacing as used in surface. + + delta = (hi-flo) / (max (ncols,nlines) -1.) * fact1 + xmin = -(real (ncols/2) * delta + real (mod (ncols+1, 2)) * delta) + ymin = -(real (nlines/2) * delta + real (mod (nlines+1, 2)) * delta) + del = 2.0 * delta + + # The perimeter is separated from the surface plot by the + # width of delta. + + x1_perim = xmin - delta + y1_perim = ymin - delta + x2_perim = xmin + (real (ncols) * delta) + y2_perim = ymin + (real (nlines) * delta) + # Set up linear arrays over full perimeter range + do i = 1, ncols + 2 + Memr[x_val+i-1] = x1_perim + (i-1) * delta + do i = 1, nlines + 2 + Memr[y_val+i-1] = y1_perim + (i-1) * delta + + # Draw and label axes and tick marks. + # It is important that frame has not been called after calling srface. + # First to draw the perimeter. Which axes get drawn depends on the + # values of angh and angv. Get angles in the range [-180, 180]. + + if (angh > 180.) + angh = angh - 360. + else if (angh < -180.) + angh = angh + 360. + + if (angv > 180.) + angv = angv - 360. + else if (angv < -180.) + angv = angv + 360. + + # Calculate positions for the axis labels + xcen = 0.5 * (x1_perim + x2_perim) + ycen = 0.5 * (y1_perim + y2_perim) + + if (angh >= 0) { + if (angv >= 0) { + # Case 1: xy rotation positive, looking down from above mid Z + + # First draw x axis + call amovkr (y2_perim, Memr[kvec], ncols + 2) + call ie_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1) + call ie_label_axis (xcen, y2_perim+del, flo, "X-AXIS", -1, -2) + call ie_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta, + flo, ncols) + junk = itoc (int (wc1), tlabel, 10) + call ie_label_axis (xmin, y2_perim+del, flo, tlabel, -1, -2) + junk = itoc (int (wc2), tlabel, 10) + call ie_label_axis (Memr[x_val+ncols], y2_perim+del, flo, + tlabel, -1, -2) + + # Now draw y axis + call amovkr (x2_perim, Memr[kvec], nlines + 2) + call ie_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1) + call ie_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1) + call ie_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1], + flo, nlines) + junk = itoc (int (wl1), tlabel, 10) + call ie_label_axis (x2_perim+del, ymin, flo, tlabel, 2, -1) + junk = itoc (int (wl2), tlabel, 10) + call ie_label_axis (x2_perim+del, Memr[y_val+nlines], flo, + tlabel, 2, -1) + } else { + # Case 2: xy rotation positive, looking up from below mid Z + # First draw x axis + call amovkr (y1_perim, Memr[kvec], ncols + 2) + call ie_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1) + call ie_label_axis (xcen, y1_perim-del, flo, "X-AXIS", -1, 2) + call ie_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, + flo, ncols) + junk = itoc (int (wc1), tlabel, 10) + call ie_label_axis (xmin, y1_perim-del, flo, tlabel, -1, 2) + junk = itoc (int (wc2), tlabel, 10) + call ie_label_axis (Memr[x_val+ncols], y1_perim-del, flo, + tlabel, -1, 2) + + # Now draw y axis + call amovkr (x1_perim, Memr[kvec], nlines + 2) + call ie_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1) + call ie_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1) + call ie_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1], + flo, nlines) + junk = itoc (int (wl1), tlabel, 10) + call ie_label_axis (x1_perim-del, ymin, flo, tlabel, 2, 1) + junk = itoc (int (wl2), tlabel, 10) + call ie_label_axis (x1_perim-del, Memr[y_val+nlines], flo, + tlabel, 2, 1) + } + } + + if (angh < 0) { + if (angv > 0) { + # Case 3: xy rotation negative, looking down from above mid Z + # (default). First draw x axis + call amovkr (y1_perim, Memr[kvec], ncols + 2) + call ie_draw_axis (Memr[x_val+1], Memr[kvec], flo, ncols + 1) + call ie_label_axis (xcen, y1_perim-del, flo, "X-AXIS", 1, 2) + call ie_draw_ticksx (Memr[x_val+1], y1_perim, y1_perim-delta, + flo, ncols) + junk = itoc (int (wc1), tlabel, 10) + call ie_label_axis (xmin, y1_perim-del, flo, tlabel, 1, 2) + junk = itoc (int (wc2), tlabel, 10) + call ie_label_axis (Memr[x_val+ncols], y1_perim-del, flo, + tlabel, 1, 2) + + # Now draw y axis + call amovkr (x2_perim, Memr[kvec], nlines + 2) + call ie_draw_axis (Memr[kvec], Memr[y_val], flo, nlines + 1) + call ie_label_axis (x2_perim+del, ycen, flo, "Y-AXIS", 2, -1) + call ie_draw_ticksy (x2_perim, x2_perim+delta, Memr[y_val+1], + flo, nlines) + junk = itoc (int (wl1), tlabel, 10) + call ie_label_axis (x2_perim+del, ymin, flo, tlabel, 2, -1) + junk = itoc (int (wl2), tlabel, 10) + call ie_label_axis (x2_perim+del, Memr[y_val+nlines], flo, + tlabel, 2, -1) + } else { + # Case 4: xy rotation negative, looking up from below mid Z + # First draw x axis + call amovkr (y2_perim, Memr[kvec], ncols + 2) + call ie_draw_axis (Memr[x_val], Memr[kvec], flo, ncols + 1) + call ie_label_axis (xcen, y2_perim+del, flo, "X-AXIS", 1, -2) + call ie_draw_ticksx (Memr[x_val+1], y2_perim, y2_perim+delta, + flo, ncols) + junk = itoc (int (wc1), tlabel, 10) + call ie_label_axis (xmin, y2_perim+del, flo, tlabel, 1, -2) + junk = itoc (int (wc2), tlabel, 10) + call ie_label_axis (Memr[x_val+ncols], y2_perim+del, flo, + tlabel, 1, -2) + + # Now draw y axis + call amovkr (x1_perim, Memr[kvec], nlines + 2) + call ie_draw_axis (Memr[kvec], Memr[y_val+1], flo, nlines + 1) + call ie_label_axis (x1_perim-del, ycen, flo, "Y-AXIS", 2, 1) + call ie_draw_ticksy (x1_perim, x1_perim-delta, Memr[y_val+1], + flo, nlines) + junk = itoc (int (wl1), tlabel, 10) + call ie_label_axis (x1_perim-del, ymin, flo, tlabel, 2, 1) + junk = itoc (int (wl2), tlabel, 10) + call ie_label_axis (x1_perim-del, Memr[y_val+nlines], flo, + tlabel, 2, 1) + } + } + + # Flush plotit buffer before returning + call plotit (0, 0, 2) + call sfree (sp) +end + + +# ?? + +procedure ie_draw_axis (xvals, yvals, zval, nvals) + +int nvals +real xvals[nvals] +real yvals[nvals] +real zval +pointer sp, xt, yt +int i +real dum + +begin + call smark (sp) + call salloc (xt, nvals, TY_REAL) + call salloc (yt, nvals, TY_REAL) + + do i = 1, nvals + call trn32s (xvals[i], yvals[i], zval, Memr[xt+i-1], Memr[yt+i-1], + dum, 1) + + call gpl (nvals, Memr[xt], Memr[yt]) + call sfree (sp) +end + + +# ?? + +procedure ie_label_axis (xval, yval, zval, sppstr, path, up) + +real xval +real yval +real zval +char sppstr[SZ_LINE] +int path +int up + +int nchars +int strlen() +% character*64 fstr + +begin + nchars = strlen (sppstr) + +% call f77pak (sppstr, fstr, 64) + call pwrzs (xval, yval, zval, fstr, nchars, CSIZE, path, up, 0) +end + + +# ?? + +procedure ie_draw_ticksx (x, y1, y2, zval, nvals) + +int nvals +real x[nvals] +real y1, y2 +real zval + +int i +real tkx[2], tky[2], dum + +begin + do i = 1, nvals { + call trn32s (x[i], y1, zval, tkx[1], tky[1], dum, 1) + call trn32s (x[i], y2, zval, tkx[2], tky[2], dum, 1) + call gpl (2, tkx[1], tky[1]) + } +end + + +# ?? + +procedure ie_draw_ticksy (x1, x2, y, zval, nvals) + +int nvals +real x1, x2 +real y[nvals] +real zval + +int i +real tkx[2], tky[2], dum + +begin + do i = 1, nvals { + call trn32s (x1, y[i], zval, tkx[1], tky[1], dum, 1) + call trn32s (x2, y[i], zval, tkx[2], tky[2], dum, 1) + call gpl (2, tkx[1], tky[1]) + } +end + + +# IE_SURF_LIMITS -- Apply the floor and ceiling constraints to the subraster. +# If either value is exactly zero, it is not applied. + +procedure ie_surf_limits (ras, m, floor, ceiling) + +real ras[m] +int m +real floor, ceiling +real val1_1 # value at ras[1] +int k +bool const_val # true if data are constant +bool bad_floor # true if no value is above floor +bool bad_ceiling # true if no value is below ceiling + +begin + const_val = true # initial values + bad_floor = true + bad_ceiling = true + val1_1 = ras[1] + + do k = 1, m + if (ras[k] != val1_1) { + const_val = false + break + } + if (!IS_INDEF(floor)) { + do k = 1, m { + if (ras[k] <= floor) + ras[k] = floor + else + bad_floor = false + } + } + if (!IS_INDEF(ceiling)) { + do k = 1, m { + if (ras[k] >= ceiling) + ras[k] = ceiling + else + bad_ceiling = false + } + } + + if (bad_floor && !IS_INDEF(floor)) + call error (1, "entire image is below (or at) specified floor") + if (bad_ceiling && !IS_INDEF(ceiling)) + call error (1, "entire image is above (or at) specified ceiling") + if (const_val) + call error (1, "all data values are the same; can't plot it") +end diff --git a/pkg/images/tv/imexamine/iestatistics.x b/pkg/images/tv/imexamine/iestatistics.x new file mode 100644 index 00000000..a3ac5f22 --- /dev/null +++ b/pkg/images/tv/imexamine/iestatistics.x @@ -0,0 +1,84 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include "imexam.h" + + +# IE_STATISTICS -- Compute and print statistics. + +procedure ie_statistics (ie, x, y) + +pointer ie # IMEXAM structure +real x, y # Aperture coordinates + +double mean, median, std +int ncstat, nlstat, x1, x2,y1, y2, npts, clgeti() +pointer sp, imname, im, data, sortdata, ie_gimage(), ie_gdata() +string label "\ +# SECTION NPIX MEAN MEDIAN STDDEV MIN MAX\n" + +begin + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + ncstat = clgeti ("ncstat") + nlstat = clgeti ("nlstat") + x1 = x - (ncstat-1) / 2 + 0.5 + x2 = x + ncstat / 2 + 0.5 + y1 = y - (nlstat-1) / 2 + 0.5 + y2 = y + nlstat / 2 + 0.5 + iferr (data = ie_gdata (im, x1, x2, y1, y2)) { + call erract (EA_WARN) + return + } + npts = (x2-x1+1) * (y2-y1+1) + + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + call salloc (sortdata, npts, TY_DOUBLE) + + call achtrd (Memr[data], Memd[sortdata], npts) + call asrtd (Memd[sortdata], Memd[sortdata], npts) + call aavgd (Memd[sortdata], npts, mean, std) + if (mod (npts, 2) == 0) + median = (Memd[sortdata+npts/2-1] + Memd[sortdata+npts/2]) / 2 + else + median = Memd[sortdata+npts/2] + + call sprintf (Memc[imname], SZ_FNAME, "[%d:%d,%d:%d]") + call pargi (x1) + call pargi (x2) + call pargi (y1) + call pargi (y2) + + if (IE_LASTKEY(ie) != 'm') + call printf (label) + + call printf ("%20s %8d %8.4g %8.4g %8.4g %8.4g %8.4g\n") + call pargstr (Memc[imname]) + call pargi (npts) + call pargd (mean) + call pargd (median) + call pargd (std) + call pargd (Memd[sortdata]) + call pargd (Memd[sortdata+npts-1]) + + if (IE_LOGFD(ie) != NULL) { + if (IE_LASTKEY(ie) != 'm') + call fprintf (IE_LOGFD(ie), label) + + call fprintf (IE_LOGFD(ie), + "%20s %8d %8.4g %8.4g %8.4g %8.4g %8.4g\n") + call pargstr (Memc[imname]) + call pargi (npts) + call pargd (mean) + call pargd (median) + call pargd (std) + call pargd (Memd[sortdata]) + call pargd (Memd[sortdata+npts-1]) + } + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/ietimexam.x b/pkg/images/tv/imexamine/ietimexam.x new file mode 100644 index 00000000..869eaa4b --- /dev/null +++ b/pkg/images/tv/imexamine/ietimexam.x @@ -0,0 +1,121 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <imhdr.h> +include "imexam.h" + + +# IE_TIMEXAM -- Extract a subraster image. +# This routine does not currently update the WCS but it does clear it. + +procedure ie_timexam (ie, x, y) + +pointer ie # IE pointer +real x, y # Center + +int i, x1, x2, y1, y2, nx, ny +pointer sp, root, extn, output +pointer im, out, data, outbuf, mw + +int clgeti(), fnextn(), iki_validextn(), strlen(), imaccess() +pointer ie_gimage(), ie_gdata(), immap(), impl2r(), mw_open() +errchk impl2r + +begin + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + call smark (sp) + call salloc (root, SZ_FNAME, TY_CHAR) + call salloc (extn, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + + # Get parameters. + call clgstr ("output", Memc[root], SZ_FNAME) + nx = clgeti ("ncoutput") + ny = clgeti ("nloutput") + + # Strip the extension. + call imgimage (Memc[root], Memc[root], SZ_FNAME) + if (Memc[root] == EOS) + call strcpy (IE_IMAGE(ie), Memc[root], SZ_FNAME) + i = fnextn (Memc[root], Memc[extn+1], SZ_FNAME) + Memc[extn] = EOS + if (i > 0) { + call iki_init() + if (iki_validextn (0, Memc[extn+1]) != 0) { + Memc[root+strlen(Memc[root])-i-1] = EOS + Memc[extn] = '.' + } + } + + do i = 1, ARB { + call sprintf (Memc[output], SZ_FNAME, "%s.%03d%s") + call pargstr (Memc[root]) + call pargi (i) + call pargstr (Memc[extn]) + if (imaccess (Memc[output], 0) == NO) + break + } + + # Set section to be extracted. + if (!IS_INDEF(x)) + IE_X1(ie) = x + if (!IS_INDEF(y)) + IE_Y1(ie) = y + + x1 = IE_X1(ie) - (nx - 1) / 2 + 0.5 + x2 = IE_X1(ie) + nx / 2 + 0.5 + y1 = IE_Y1(ie) - (ny - 1) / 2 + 0.5 + y2 = IE_Y1(ie) + ny / 2 + 0.5 + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + + # Set output. + iferr (out = immap (Memc[output], NEW_COPY, im)) { + call erract (EA_WARN) + return + } + IM_NDIM(out) = 2 + IM_LEN(out,1) = nx + IM_LEN(out,2) = ny + + # Extract the section. + iferr { + do i = y1, y2 { + data = ie_gdata (im, x1, x2, i, i) + outbuf = impl2r (out, i-y1+1) + call amovr (Memr[data], Memr[outbuf], nx) + } + mw = mw_open (NULL, 2) + call mw_saveim (mw, out) + call imunmap (out) + } then { + call imunmap (out) + iferr (call imdelete (Memc[output])) + ; + call sfree (sp) + call erract (EA_WARN) + return + } + + call printf ("%s[%d:%d,%d:%d] -> %s\n") + call pargstr (IE_IMAGE(ie)) + call pargi (x1) + call pargi (x2) + call pargi (y1) + call pargi (y2) + call pargstr (Memc[output]) + if (IE_LOGFD(ie) != NULL) { + call fprintf (IE_LOGFD(ie), "%s[%d:%d,%d:%d] -> %s\n") + call pargstr (IE_IMAGE(ie)) + call pargi (x1) + call pargi (x2) + call pargi (y1) + call pargi (y2) + } + + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/ievimexam.x b/pkg/images/tv/imexamine/ievimexam.x new file mode 100644 index 00000000..a75ac2bc --- /dev/null +++ b/pkg/images/tv/imexamine/ievimexam.x @@ -0,0 +1,582 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <gset.h> +include <mach.h> +include <math.h> +include <imhdr.h> +include <imset.h> +include <math/iminterp.h> +include "imexam.h" + +define BTYPES "|constant|nearest|reflect|wrap|project|" +define SZ_BTYPE 8 # Length of boundary type string +define NLINES 16 # Number of image lines in the buffer + + +# IE_VIMEXAM -- Plot the vector of image data between two pixels. +# There are two types of plot selected by the key argument. The +# second cursor position is passed in the IMEXAM data structure. +# The first position is either the middle of the vector or the starting +# point. + +procedure ie_vimexam (gp, mode, ie, x, y, key) + +pointer gp # GIO pointer +int mode # Graph mode +pointer ie # IMEXAM pointer +real x, y # Starting or center coordinate +int key # 'u' centered vector, 'v' two endpoint vector + +int btype, nxvals, nyvals, nzvals, width +pointer sp, title, boundary, im, x_vec, y_vec, pp +real x1, y1, x2, y2, zmin, zmax, bconstant + +bool fp_equalr() +int clgpseti(), clgwrd(), clopset() +real clgpsetr() +pointer ie_gimage() +errchk malloc + +begin + iferr (im = ie_gimage (ie, NO)) { + call erract (EA_WARN) + return + } + + call smark (sp) + call salloc (title, IE_SZTITLE, TY_CHAR) + call salloc (boundary, SZ_BTYPE, TY_CHAR) + + # Get boundary extension parameters. + if (IE_PP(ie) != NULL) + call clcpset (IE_PP(ie)) + IE_PP(ie) = clopset ("vimexam") + pp = IE_PP(ie) + btype = clgwrd ("vimexam.boundary", Memc[boundary], SZ_BTYPE, BTYPES) + bconstant = clgpsetr (pp, "constant") + + nxvals = IM_LEN(im,1) + nyvals = IM_LEN(im,2) + + if (!IS_INDEF (x)) + IE_X1(ie) = x + if (!IS_INDEF(y)) + IE_Y1(ie) = y + + x1 = IE_X1(ie) + x2 = IE_X2(ie) + y1 = IE_Y1(ie) + y2 = IE_Y2(ie) + width = clgpseti (pp, "naverage") + + # Check the boundary and compute the length of the output vector. + x1 = max (1.0, min (x1, real (nxvals))) + x2 = min (real(nxvals), max (1.0, x2)) + y1 = max (1.0, min (y1, real (nyvals))) + y2 = min (real(nyvals), max (1.0, y2)) + nzvals = int (sqrt ((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1))) + 1 + + # Check for cases which should be handled by pcols or prows. + call malloc (x_vec, nzvals, TY_REAL) + call malloc (y_vec, nzvals, TY_REAL) + if (fp_equalr (x1, x2)) + call ie_get_col (im, x1, y1, x2, y2, nzvals, width, btype, + bconstant, Memr[x_vec], Memr[y_vec], zmin, zmax) + else if (fp_equalr (y1, y2)) + call ie_get_row (im, x1, y1, x2, y2, nzvals, width, btype, + bconstant, Memr[x_vec], Memr[y_vec], zmin, zmax) + else + call ie_get_vector (im, x1, y1, x2, y2, nzvals, width, btype, + bconstant, Memr[x_vec], Memr[y_vec], zmin, zmax) + + # Convert endpoint plot coordinates to centered coordinates. + if (key == 'u') { + zmin = (IE_X1(ie) + IE_X2(ie)) / 2 + zmax = (IE_Y1(ie) + IE_Y2(ie)) / 2 + zmin = sqrt ((zmin-x1)**2 + (zmax-y1)**2) + call asubkr (Memr[x_vec], zmin, Memr[x_vec], nzvals) + } + + call sprintf (Memc[title], IE_SZTITLE, + "%s: Vector %.1f,%.1f to %.1f,%.1f naverage: %d\n%s") + call pargstr (IE_IMNAME(ie)) + call pargr (x1) + call pargr (y1) + call pargr (x2) + call pargr (y2) + call pargi (width) + call pargstr (IM_TITLE(im)) + + call ie_graph (gp, mode, pp, Memc[title], Memr[x_vec], Memr[y_vec], + nzvals, "", "") + + # Finish up + call mfree (x_vec, TY_REAL) + call mfree (y_vec, TY_REAL) + call sfree (sp) +end + + +# IE_GET_VECTOR -- Average a strip perpendicular to a given vector and return +# vectors of point number and average pixel value. Also returned is the min +# and max value in the data vector. + +procedure ie_get_vector (im, x1, y1, x2, y2, nvals, width, btype, + bconstant, x_vector, y_vector, zmin, zmax) + +pointer im # pointer to image header +real x1, y1 # starting pixel of vector +real x2, y2 # ending pixel of pixel +real bconstant # Boundary extension constant +int btype # Boundary extension type +int nvals # number of samples along the vector +int width # width of strip to average over +real x_vector[ARB] # Pixel numbers +real y_vector[ARB] # Average pixel values (returned) +real zmin, zmax # min, max of data vector + +double dx, dy, dpx, dpy, ratio, xoff, yoff, noff, xv, yv +int i, j, k, nedge, col1, col2, line1, line2 +int colb, colc, line, linea, lineb, linec +pointer sp, oxs, oys, xs, ys, yvals, msi, buf +real sum , lim1, lim2, lim3, lim4 +pointer imgs2r() +errchk msiinit + +begin + call smark (sp) + call salloc (oxs, width, TY_REAL) + call salloc (oys, width, TY_REAL) + call salloc (xs, width, TY_REAL) + call salloc (ys, width, TY_REAL) + call salloc (yvals, width, TY_REAL) + + # Determine sampling perpendicular to vector. + dx = (x2 - x1) / (nvals - 1) + dy = (y2 - y1) / (nvals - 1) + if (x1 < x2) { + dpx = -dy + dpy = dx + } else { + dpx = dy + dpy = -dx + } + + # Compute offset from the nominal vector to the first sample point. + ratio = dx / dy + nedge = width + 1 + noff = (real (width) - 1.0) / 2.0 + xoff = noff * dpx + yoff = noff * dpy + + # Initialize the interpolator and the image data buffer. + call msiinit (msi, II_BILINEAR) + buf = NULL + + # Set the boundary. + col1 = int (min (x1, x2)) - nedge + col2 = nint (max (x1, x2)) + nedge + line1 = int (min (y1, y2)) - nedge + line2 = nint (max (y2, y1)) + nedge + call ie_setboundary (im, col1, col2, line1, line2, btype, bconstant) + + # Initialize. + xv = x1 - xoff + yv = y1 - yoff + do j = 1, width { + Memr[oxs+j-1] = double (j - 1) * dpx + Memr[oys+j-1] = double (j - 1) * dpy + } + + # Loop over the output image lines. + do i = 1, nvals { + x_vector[i] = real (i) + line = yv + + # Get the input image data and fit an interpolator to the data. + # The input data is buffered in a section of size NLINES + 2 * + # NEDGE. + + if (dy >= 0.0 && (buf == NULL || line > (linea))) { + linea = min (line2, line + NLINES - 1) + lineb = max (line1, line - nedge) + linec = min (line2, linea + nedge) + lim1 = xv + lim2 = lim1 + double (width - 1) * dpx + lim3 = xv + double (linea - line + 1) * ratio + lim4 = lim3 + double (width - 1) * dpx + colb = max (col1, int (min (lim1, lim2, lim3, lim4)) - 1) + colc = min (col2, nint (max (lim1, lim2, lim3, lim4)) + 1) + buf = imgs2r (im, colb, colc, lineb, linec) + call msifit (msi, Memr[buf], colc - colb + 1, linec - lineb + + 1, colc - colb + 1) + + } else if (dy < 0.0 && (buf == NULL || line < linea)) { + linea = max (line1, line - NLINES + 1) + lineb = max (line1, linea - nedge) + linec = min (line2, line + nedge) + lim1 = xv + lim2 = lim1 + double (width - 1) * dpx + lim3 = xv + double (linea - line - 1) * ratio + lim4 = lim3 + double (width - 1) * dpx + colb = max (col1, int (min (lim1, lim2, lim3, lim4)) - 1) + colc = min (col2, nint (max (lim1, lim2, lim3, lim4)) + 1) + buf = imgs2r (im, colb, colc, lineb, linec) + call msifit (msi, Memr[buf], colc - colb + 1, linec - lineb + + 1, colc - colb + 1) + } + + # Evaluate the interpolant. + call aaddkr (Memr[oxs], real (xv - colb + 1), Memr[xs], width) + call aaddkr (Memr[oys], real (yv - lineb + 1), Memr[ys], width) + call msivector (msi, Memr[xs], Memr[ys], Memr[yvals], width) + + if (width == 1) + y_vector[i] = Memr[yvals] + else { + sum = 0.0 + do k = 1, width + sum = sum + Memr[yvals+k-1] + y_vector[i] = sum / width + } + + xv = xv + dx + yv = yv + dy + } + + # Compute min and max values. + call alimr (y_vector, nvals, zmin, zmax) + + # Free memory . + call msifree (msi) + call sfree (sp) +end + + +# IE_GET_COL -- Average a strip perpendicular to a column vector and return +# vectors of point number and average pixel value. Also returned is the min +# and max value in the data vector. + +procedure ie_get_col (im, x1, y1, x2, y2, nvals, width, btype, + bconstant, x_vector, y_vector, zmin, zmax) + +pointer im # pointer to image header +real x1, y1 # starting pixel of vector +real x2, y2 # ending pixel of pixel +int nvals # number of samples along the vector +int width # width of strip to average over +int btype # Boundary extension type +real bconstant # Boundary extension constant +real x_vector[ARB] # Pixel numbers +real y_vector[ARB] # Average pixel values (returned) +real zmin, zmax # min, max of data vector + +real sum +int line, linea, lineb, linec +pointer sp, xs, ys, msi, yvals, buf +double dx, dy, xoff, noff, xv, yv +int i, j, k, nedge, col1, col2, line1, line2 +pointer imgs2r() +errchk msiinit + +begin + call smark (sp) + call salloc (xs, width, TY_REAL) + call salloc (ys, width, TY_REAL) + call salloc (yvals, width, TY_REAL) + + # Initialize the interpolator and the image data buffer. + call msiinit (msi, II_BILINEAR) + buf = NULL + + # Set the boundary. + nedge = max (2, width / 2 + 1) + col1 = int (x1) - nedge + col2 = nint (x1) + nedge + line1 = int (min (y1, y2)) - nedge + line2 = nint (max (y1, y2)) + nedge + call ie_setboundary (im, col1, col2, line1, line2, btype, bconstant) + + # Determine sampling perpendicular to vector. + dx = 1.0d0 + if (nvals == 1) + dy = 0.0d0 + else + dy = (y2 - y1) / (nvals - 1) + + # Compute offset from the nominal vector to the first sample point. + noff = (real (width) - 1.0) / 2.0 + xoff = noff * dx + xv = x1 - xoff + do j = 1, width + Memr[xs+j-1] = xv + double (j - col1) + yv = y1 + + # Loop over the output image lines. + do i = 1, nvals { + x_vector[i] = real (i) + line = yv + + # Get the input image data and fit an interpolator to the data. + # The input data is buffered in a section of size NLINES + 2 * + # NEDGE. + + if (dy >= 0.0 && (buf == NULL || line > (linea))) { + linea = min (line2, line + NLINES - 1) + lineb = max (line1, line - nedge) + linec = min (line2, linea + nedge) + buf = imgs2r (im, col1, col2, lineb, linec) + call msifit (msi, Memr[buf], col2 - col1 + 1, linec - lineb + + 1, col2 - col1 + 1) + } else if (dy < 0.0 && (buf == NULL || line < linea)) { + linea = max (line1, line - NLINES + 1) + lineb = max (line1, linea - nedge) + linec = min (line2, line + nedge) + buf = imgs2r (im, col1, col2, lineb, linec) + call msifit (msi, Memr[buf], col2 - col1 + 1, linec - lineb + + 1, col2 - col1 + 1) + } + + # Evaluate the interpolant. + call amovkr (real (yv - lineb + 1), Memr[ys], width) + call msivector (msi, Memr[xs], Memr[ys], Memr[yvals], width) + + if (width == 1) + y_vector[i] = Memr[yvals] + else { + sum = 0.0 + do k = 1, width + sum = sum + Memr[yvals+k-1] + y_vector[i] = sum / width + } + + yv = yv + dy + } + + # Compute min and max values. + call alimr (y_vector, nvals, zmin, zmax) + + # Free memory . + call msifree (msi) + call sfree (sp) +end + + +# IE_GET_ROW -- Average a strip parallel to a row vector and return +# vectors of point number and average pixel value. Also returned is the min +# and max value in the data vector. + +procedure ie_get_row (im, x1, y1, x2, y2, nvals, width, btype, bconstant, + x_vector, y_vector, zmin, zmax) + +pointer im # pointer to image header +real x1, y1 # starting pixel of vector +real x2, y2 # ending pixel of pixel +int nvals # number of samples along the vector +int width # width of strip to average over +int btype # Boundary extension type +real bconstant # Boundary extension constant +real x_vector[ARB] # Pixel numbers +real y_vector[ARB] # Average pixel values (returned) +real zmin, zmax # min, max of data vector + +double dx, dy, yoff, noff, xv, yv +int i, j, nedge, col1, col2, line1, line2 +int line, linea, lineb, linec +pointer sp, oys, xs, ys, yvals, msi, buf +errchk imgs2r, msifit, msiinit +pointer imgs2r() + +begin + call smark (sp) + call salloc (oys, width, TY_REAL) + call salloc (xs, nvals, TY_REAL) + call salloc (ys, nvals, TY_REAL) + call salloc (yvals, nvals, TY_REAL) + + # Initialize the interpolator and the image data buffer. + call msiinit (msi, II_BILINEAR) + buf = NULL + + # Set the boundary. + nedge = max (2, width / 2 + 1) + col1 = int (min (x1, x2)) - nedge + col2 = nint (max (x1, x2)) + nedge + line1 = int (y1) - nedge + line2 = nint (y1) + nedge + call ie_setboundary (im, col1, col2, line1, line2, btype, bconstant) + + # Determine sampling perpendicular to vector. + if (nvals == 1) + dx = 0.0d0 + else + dx = (x2 - x1) / (nvals - 1) + dy = 1.0 + + # Compute offset from the nominal vector to the first sample point. + noff = (real (width) - 1.0) / 2.0 + xv = x1 - col1 + 1 + do i = 1, nvals { + Memr[xs+i-1] = xv + xv = xv + dx + } + yoff = noff * dy + yv = y1 - yoff + do j = 1, width + Memr[oys+j-1] = yv + double (j - 1) + + # Clear the accululator. + call aclrr (y_vector, nvals) + + # Loop over the output image lines. + do i = 1, width { + line = yv + + # Get the input image data and fit an interpolator to the data. + # The input data is buffered in a section of size NLINES + 2 * + # NEDGE. + + if (dy >= 0.0 && (buf == NULL || line > (linea))) { + linea = min (line2, line + NLINES - 1) + lineb = max (line1, line - nedge) + linec = min (line2, linea + nedge) + buf = imgs2r (im, col1, col2, lineb, linec) + if (buf == NULL) + call error (0, "Error reading input image.") + call msifit (msi, Memr[buf], col2 - col1 + 1, linec - lineb + + 1, col2 - col1 + 1) + } else if (dy < 0.0 && (buf == NULL || line < linea)) { + linea = max (line1, line - NLINES + 1) + lineb = max (line1, linea - nedge) + linec = min (line2, line + nedge) + buf = imgs2r (im, col1, col2, lineb, linec) + if (buf == NULL) + call error (0, "Error reading input image.") + call msifit (msi, Memr[buf], col2 - col1 + 1, linec - lineb + + 1, col2 - col1 + 1) + } + + # Evaluate the interpolant. + call amovkr (real (Memr[oys+i-1] - lineb + 1), Memr[ys], nvals) + call msivector (msi, Memr[xs], Memr[ys], Memr[yvals], nvals) + + if (width == 1) + call amovr (Memr[yvals], y_vector, nvals) + else + call aaddr (Memr[yvals], y_vector, y_vector, nvals) + + yv = yv + dy + } + + # Compute the x and y vectors. + do i = 1, nvals + x_vector[i] = real (i) + if (width > 1) + call adivkr (y_vector, real (width), y_vector, nvals) + + # Compute min and max values. + call alimr (y_vector, nvals, zmin, zmax) + + # Free memory . + call msifree (msi) + call sfree (sp) +end + + +# IE_SETBOUNDARY -- Set boundary extension. + +procedure ie_setboundary (im, col1, col2, line1, line2, btype, bconstant) + +pointer im # IMIO pointer +int col1, col2 # Range of columns +int line1, line2 # Range of lines +int btype # Boundary extension type +real bconstant # Constant for constant boundary extension + +int btypes[5] +int nbndrypix +data btypes /BT_CONSTANT, BT_NEAREST, BT_REFLECT, BT_WRAP, BT_PROJECT/ + +begin + nbndrypix = 0 + nbndrypix = max (nbndrypix, 1 - col1) + nbndrypix = max (nbndrypix, col2 - IM_LEN(im, 1)) + nbndrypix = max (nbndrypix, 1 - line1) + nbndrypix = max (nbndrypix, line2 - IM_LEN(im, 2)) + + call imseti (im, IM_TYBNDRY, btypes[btype]) + call imseti (im, IM_NBNDRYPIX, nbndrypix + 1) + if (btypes[btype] == BT_CONSTANT) + call imsetr (im, IM_BNDRYPIXVAL, bconstant) +end + + +# IE_BUFL2R -- Maintain buffer of image lines. A new buffer is created when +# the buffer pointer is null or if the number of lines requested is changed. +# The minimum number of image reads is used. + +procedure ie_bufl2r (im, col1, col2, line1, line2, buf) + +pointer im # Image pointer +int col1 # First image column of buffer +int col2 # Last image column of buffer +int line1 # First image line of buffer +int line2 # Last image line of buffer +pointer buf # Buffer + +pointer buf1, buf2 +int i, ncols, nlines, nclast, llast1, llast2, nllast +errchk malloc, realloc, imgs2r +pointer imgs2r() + +begin + ncols = col2 - col1 + 1 + nlines = line2 - line1 + 1 + + # If the buffer pointer is undefined then allocate memory for the + # buffer. If the number of columns or lines requested changes + # reallocate the buffer. Initialize the last line values to force + # a full buffer image read. + + if (buf == NULL) { + call malloc (buf, ncols * nlines, TY_REAL) + llast1 = line1 - nlines + llast2 = line2 - nlines + } else if ((nlines != nllast) || (ncols != nclast)) { + call realloc (buf, ncols * nlines, TY_REAL) + llast1 = line1 - nlines + llast2 = line2 - nlines + } + + # Read only the image lines with are different from the last buffer. + + if (line1 < llast1) { + do i = line2, line1, -1 { + if (i > llast1) + buf1 = buf + (i - llast1) * ncols + else + buf1 = imgs2r (im, col1, col2, i, i) + + buf2 = buf + (i - line1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } else if (line2 > llast2) { + do i = line1, line2 { + if (i < llast2) + buf1 = buf + (i - llast1) * ncols + else + buf1 = imgs2r (im, col1, col2, i, i) + + buf2 = buf + (i - line1) * ncols + call amovr (Memr[buf1], Memr[buf2], ncols) + } + } + + # Save the buffer parameters. + + llast1 = line1 + llast2 = line2 + nclast = ncols + nllast = nlines +end diff --git a/pkg/images/tv/imexamine/imexam.h b/pkg/images/tv/imexamine/imexam.h new file mode 100644 index 00000000..f1fe00d8 --- /dev/null +++ b/pkg/images/tv/imexamine/imexam.h @@ -0,0 +1,55 @@ +# IMEXAM.H -- IMEXAMINE global definitions. + +define MAX_FRAMES 16 # max display frames + +# IMEXAMINE data structure. + +define IE_LEN 370 # length of IE structure +define IE_SZFNAME 99 # length of file name +define IE_SZFORMAT 9 # length of format strings +define IE_SZTITLE 512 # length of multiline title + +define IE_IM Memi[$1] # IMIO pointer +define IE_MW Memi[$1+1] # MWCS pointer +define IE_CTLW Memi[$1+2] # CT-MWCS pointer (L -> W) +define IE_CTWL Memi[$1+3] # CT-MWCS pointer (W -> L) +define IE_DS Memi[$1+4] # display frame pointer +define IE_GP Memi[$1+5] # GIO pointer +define IE_PP Memi[$1+6] # pset pointer +define IE_LIST Memi[$1+7] # image list +define IE_LISTLEN Memi[$1+8] # number of images in list +define IE_USEDISPLAY Memi[$1+9] # use image display? +define IE_INDEX Memi[$1+10] # image index +define IE_DFRAME Memi[$1+11] # frame used to display images +define IE_MAPFRAME Memi[$1+12] # mapped display frame +define IE_NEWFRAME Memi[$1+13] # new (current) display frame +define IE_NFRAMES Memi[$1+14] # number of image frames +define IE_ALLFRAMES Memi[$1+15] # use all frames for display? +define IE_LOGFD Memi[$1+16] # log file descriptor +define IE_MAGZERO Memr[P2R($1+17)] # magnitude zero point +define IE_XORIGIN Memr[P2R($1+18)] # X origin +define IE_YORIGIN Memr[P2R($1+19)] # Y origin +define IE_GTYPE Memi[$1+20] # current graph type +define IE_X1 Memr[P2R($1+21)] # current graph x1 +define IE_X2 Memr[P2R($1+22)] # current graph x2 +define IE_Y1 Memr[P2R($1+23)] # current graph y1 +define IE_Y2 Memr[P2R($1+24)] # current graph y2 +define IE_IX1 Memi[$1+25] # image section coordinate +define IE_IX2 Memi[$1+26] # image section coordinate +define IE_IY1 Memi[$1+27] # image section coordinate +define IE_IY2 Memi[$1+28] # image section coordinate +define IE_P1 Memi[$1+29] # Physical axis for logical x +define IE_P2 Memi[$1+30] # Physical axis for logical y +define IE_IN Memr[P2R($1+31)+$2-1] # Input coordinate vector +define IE_OUT Memr[P2R($1+38)+$2-1] # Output coordinate vector +define IE_WCSDIM Memi[$1+45] # WCS dimension +define IE_LASTKEY Memi[$1+46] # last type of keyed output + # (available) +define IE_IMAGE Memc[P2C($1+50)] # full image name +define IE_IMNAME Memc[P2C($1+100)] # short image name for labels +define IE_LOGFILE Memc[P2C($1+150)] # logfile name +define IE_WCSNAME Memc[P2C($1+200)] # WCS name +define IE_XLABEL Memc[P2C($1+250)] # WCS label +define IE_YLABEL Memc[P2C($1+300)] # WCS label +define IE_XFORMAT Memc[P2C($1+350)] # WCS format +define IE_YFORMAT Memc[P2C($1+360)] # WCS format diff --git a/pkg/images/tv/imexamine/imexamine.par b/pkg/images/tv/imexamine/imexamine.par new file mode 100644 index 00000000..fc409b45 --- /dev/null +++ b/pkg/images/tv/imexamine/imexamine.par @@ -0,0 +1,22 @@ +input,s,a,,,,images to be examined +output,s,h,"",,,output root image name +ncoutput,i,h,101,1,,Number of columns in image output +nloutput,i,h,101,1,,Number of lines in image output +frame,i,q,1,1,,display frame +image,s,q,,,,image name +logfile,s,h,"",,,logfile +keeplog,b,h,no,,,log output results +defkey,s,h,"a",,,default key for cursor list input +autoredraw,b,h,yes,,,automatically redraw graph +allframes,b,h,yes,,,use all frames for displaying new images +nframes,i,h,0,,,number of display frames (0 to autosense) +ncstat,i,h,5,1,,number of columns for statistics +nlstat,i,h,5,1,,number of lines for statistics +graphcur,*gcur,h,"",,,graphics cursor input +imagecur,*imcur,h,"",,,image display cursor input +wcs,s,h,"logical",,,Coordinate system +xformat,s,h,"",,,X axis coordinate format +yformat,s,h,"",,,Y axis coordinate format +graphics,s,h,"stdgraph",,,graphics device +display,s,h,"display(image='$1',frame=$2)",,,display command template +use_display,b,h,yes,,,enable direct display interaction diff --git a/pkg/images/tv/imexamine/mkpkg b/pkg/images/tv/imexamine/mkpkg new file mode 100644 index 00000000..38c3fef7 --- /dev/null +++ b/pkg/images/tv/imexamine/mkpkg @@ -0,0 +1,48 @@ +# IMEXAMINE + +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +standalone: + $set LIBS1 = "-lds -liminterp -lncar -lgks -lxtools" + $set LIBS2 = "-lgsurfit -lnlfit -lcurfit -lllsq" + $update libpkg.a + $omake x_imexam.x + $link x_imexam.o libpkg.a $(LIBS1) $(LIBS2) -o xx_imexam.e + ; + +libpkg.a: + iecimexam.x imexam.h <error.h> <imhdr.h> + iecolon.x imexam.h <error.h> <imhdr.h> + iedisplay.x <error.h> + ieeimexam.x imexam.h <config.h> <error.h> <fset.h> <gset.h>\ + <imhdr.h> <mach.h> <xwhen.h> + iegcur.x imexam.h <imhdr.h> <ctype.h> <mach.h> + iegdata.x <imhdr.h> + iegimage.x imexam.h <error.h> <imhdr.h> + iegnfr.x imexam.h <imhdr.h> + iegraph.x imexam.h <gset.h> + iehimexam.x imexam.h <error.h> <imhdr.h> + ieimname.x + iejimexam.x imexam.h <error.h> <imhdr.h> <gset.h> <mach.h> + ielimexam.x imexam.h <error.h> <imhdr.h> + iemw.x imexam.h <imhdr.h> <mwset.h> + ieopenlog.x imexam.h <error.h> <imhdr.h> + iepos.x imexam.h <error.h> <math.h> + ieprint.x imexam.h <error.h> + ieqrimexam.x imexam.h <error.h> <imhdr.h> <gset.h> <math.h>\ + <math/gsurfit.h> <math/nlfit.h> + ierimexam.x imexam.h <error.h> <gset.h> <imhdr.h> <math.h>\ + <math/gsurfit.h> <math/nlfit.h> + iesimexam.x imexam.h <error.h> <gset.h> <imhdr.h> <mach.h> + iestatistics.x imexam.h <error.h> + ietimexam.x imexam.h <error.h> <imhdr.h> + ievimexam.x imexam.h <error.h> <gset.h> <imhdr.h> <mach.h>\ + <imset.h> <math.h> <math/iminterp.h> + stfmeasure.x starfocus.h <error.h> <imhdr.h> <imset.h> <math/nlfit.h> + stfprofile.x starfocus.h <imhdr.h> <mach.h>\ + <math.h> <math/nlfit.h> <math/iminterp.h> + t_imexam.x imexam.h <error.h> <gset.h> <imhdr.h> + ; diff --git a/pkg/images/tv/imexamine/starfocus.h b/pkg/images/tv/imexamine/starfocus.h new file mode 100644 index 00000000..cf397e50 --- /dev/null +++ b/pkg/images/tv/imexamine/starfocus.h @@ -0,0 +1,140 @@ +# STARFOCUS + +# Types of coordinates +define SF_TYPES "|center|mark1|markall|" +define SF_CENTER 1 # Star at center of image +define SF_MARK1 2 # Mark stars in first image +define SF_MARKALL 3 # Mark stars in all images + +# Task type +define STARFOCUS 1 +define PSFMEASURE 2 + +# Radius types +define SF_WTYPES "|Radius|FWHM|GFWHM|MFWHM|" + +define SF_RMIN 16 # Minimum centering search radius +define MAX_FRAMES 8 # Maximum number of display frames + +# Data structures for STARFOCUS + +define NBNDRYPIX 0 # Number of boundary pixels +define TYBNDRY BT_REFLECT # Type of boundary extension +define SAMPLE .2 # Subpixel sampling size +define SF_SZFNAME 79 # Length of file names +define SF_SZWTYPE 7 # Length of width type string + +# Main data structure +define SF 40 +define SF_TASK Memi[$1] # Task type +define SF_WTYPE Memc[P2C($1+1)] # Width type string +define SF_WCODE Memi[$1+5] # Width code +define SF_BETA Memr[P2R($1+6)] # Moffat beta +define SF_SCALE Memr[P2R($1+7)] # Pixel scale +define SF_LEVEL Memr[P2R($1+8)] # Profile measurement level +define SF_RADIUS Memr[P2R($1+9)] # Profile radius +define SF_SBUF Memr[P2R($1+10)] # Sky region buffer +define SF_SWIDTH Memr[P2R($1+11)] # Sky region width +define SF_SAT Memr[P2R($1+12)] # Saturation +define SF_NIT Memi[$1+13] # Number of iterations for radius +define SF_OVRPLT Memi[$1+14] # Overplot the best profile? +define SF_NCOLS Memi[$1+15] # Number of image columns +define SF_NLINES Memi[$1+16] # Number of image lines +define SF_XF Memr[P2R($1+17)] # X field center +define SF_YF Memr[P2R($1+18)] # Y field center +define SF_GP Memi[$1+19] # GIO pointer +define SF_F Memr[P2R($1+20)] # Best focus +define SF_W Memr[P2R($1+21)] # Width at best focus +define SF_M Memr[P2R($1+22)] # Brightest star magnitude +define SF_XP1 Memr[P2R($1+23)] # First derivative point to plot +define SF_XP2 Memr[P2R($1+24)] # Last derivative point to plot +define SF_YP1 Memr[P2R($1+25)] # Minimum of derivative profile +define SF_YP2 Memr[P2R($1+26)] # Maximum of derivative profile +define SF_N Memi[$1+27] # Number of points not deleted +define SF_NSFD Memi[$1+28] # Number of data points +define SF_SFDS Memi[$1+29] # Pointer to data structures +define SF_NS Memi[$1+30] # Number of stars not deleted +define SF_NSTARS Memi[$1+31] # Number of stars +define SF_STARS Memi[$1+32] # Pointer to star groups +define SF_NF Memi[$1+33] # Number of focuses not deleted +define SF_NFOCUS Memi[$1+34] # Number of different focus values +define SF_FOCUS Memi[$1+35] # Pointer to focus groups +define SF_NI Memi[$1+36] # Number of images not deleted +define SF_NIMAGES Memi[$1+37] # Number of images +define SF_IMAGES Memi[$1+38] # Pointer to image groups +define SF_BEST Memi[$1+39] # Pointer to best focus star + +define SF_SFD Memi[SF_SFDS($1)+$2-1] +define SF_SFS Memi[SF_STARS($1)+$2-1] +define SF_SFF Memi[SF_FOCUS($1)+$2-1] +define SF_SFI Memi[SF_IMAGES($1)+$2-1] + +# Basic data structure. +define SFD 94 +define SFD_IMAGE Memc[P2C($1)] # Image name +define SFD_DATA Memi[$1+40] # Pointer to real image raster +define SFD_RADIUS Memr[P2R($1+41)] # Profile radius +define SFD_NP Memi[$1+42] # Number of profile points +define SFD_NPMAX Memi[$1+43] # Maximum number of profile points +define SFD_X1 Memi[$1+44] # Image raster limits +define SFD_X2 Memi[$1+45] +define SFD_Y1 Memi[$1+46] +define SFD_Y2 Memi[$1+47] +define SFD_ID Memi[$1+48] # Star ID +define SFD_X Memr[P2R($1+49)] # Star X position +define SFD_Y Memr[P2R($1+50)] # Star Y position +define SFD_F Memr[P2R($1+51)] # Focus +define SFD_W Memr[P2R($1+52)] # Width to use +define SFD_M Memr[P2R($1+53)] # Magnitude +define SFD_E Memr[P2R($1+54)] # Ellipticity +define SFD_PA Memr[P2R($1+55)] # Position angle +define SFD_R Memr[P2R($1+56)] # Radius at given level +define SFD_DFWHM Memr[P2R($1+57)] # Direct FWHM +define SFD_GFWHM Memr[P2R($1+58)] # Gaussian FWHM +define SFD_MFWHM Memr[P2R($1+59)] # Moffat FWHM +define SFD_ASI1 Memi[$1+60] # Pointer to enclosed flux profile +define SFD_ASI2 Memi[$1+61] # Pointer to derivative profile +define SFD_YP1 Memr[P2R($1+62)] # Minimum of derivative profile +define SFD_YP2 Memr[P2R($1+63)] # Maximum of derivative profile +define SFD_FWHM Memr[P2R($1+$2+63)] # FWHM vs level=0.5*i (i=1-19) +define SFD_BKGD Memr[P2R($1+83)] # Background value +define SFD_BKGD1 Memr[P2R($1+84)] # Original background value +define SFD_MISO Memr[P2R($1+85)] # Moment isophote +define SFD_SIGMA Memr[P2R($1+86)] # Moffat alpha +define SFD_ALPHA Memr[P2R($1+87)] # Moffat alpha +define SFD_BETA Memr[P2R($1+88)] # Moffat beta +define SFD_STATUS Memi[$1+89] # Status +define SFD_NSAT Memi[$1+90] # Number of saturated pixels +define SFD_SFS Memi[$1+91] # Pointer to star group +define SFD_SFF Memi[$1+92] # Pointer to focus group +define SFD_SFI Memi[$1+93] # Pointer to image group + + +# Structure grouping data by star. +define SFS ($1+7) +define SFS_ID Memi[$1] # Star ID +define SFS_F Memr[P2R($1+1)] # Best focus +define SFS_W Memr[P2R($1+2)] # Best width +define SFS_M Memr[P2R($1+3)] # Average magnitude +define SFS_N Memi[$1+4] # Number of points used +define SFS_NF Memi[$1+5] # Number of focuses +define SFS_NSFD Memi[$1+6] # Number of data points +define SFS_SFD Memi[$1+$2+6] # Array of data structures + + +# Structure grouping stars by focus values. +define SFF ($1+5) +define SFF_F Memr[P2R($1)] # Focus +define SFF_W Memr[P2R($1+1)] # Average width +define SFF_N Memi[$1+2] # Number in average +define SFF_NI Memi[$1+3] # Number of images +define SFF_NSFD Memi[$1+4] # Number of data points +define SFF_SFD Memi[$1+$2+4] # Array of data structures + + +# Structure grouping stars by image. +define SFI ($1+42) +define SFI_IMAGE Memc[P2C($1)] # Image +define SFI_N Memi[$1+40] # Number in imagE +define SFI_NSFD Memi[$1+41] # Number of data points +define SFI_SFD Memi[$1+$2+41] # Array of data structures diff --git a/pkg/images/tv/imexamine/stfmeasure.x b/pkg/images/tv/imexamine/stfmeasure.x new file mode 100644 index 00000000..7390bf1c --- /dev/null +++ b/pkg/images/tv/imexamine/stfmeasure.x @@ -0,0 +1,147 @@ +include <error.h> +include <imhdr.h> +include <imset.h> +include <math/nlfit.h> +include "starfocus.h" + + +# STF_MEASURE -- PSF measuring routine. +# This is a stand-alone routine that can be called to return the FWHM. +# It is a greatly abbreviated version of starfocus. + +procedure stf_measure (im, xc, yc, beta, level, radius, nit, + sbuffer, swidth, saturation, gp, logfd, + bkg, renclosed, dfwhm, gfwhm, mfwhm) + +pointer im #I Image pointer +real xc #I Initial X center +real yc #I Initial Y center +real beta #I Moffat beta +real level #I Measurement level +real radius #U Profile radius +int nit #I Number of iterations on radius +real sbuffer #I Sky buffer (pixels) +real swidth #I Sky width (pixels) +real saturation #I Saturation +pointer gp #I Graphics output if not NULL +int logfd #I Log output if not NULL +real bkg #O Background used +real renclosed #O Enclosed flux radius +real dfwhm #O Direct FWHM +real gfwhm #O Gaussian FWHM +real mfwhm #O Moffat FWHM + +int i +bool ignore_sat +pointer sp, str, sf, sfd, sfds + +int strdic() +real stf_r2i() +errchk stf_find, stf_bkgd, stf_profile, stf_widths, stf_fwhms, stf_organize + +begin + call smark (sp) + call salloc (str, SZ_FNAME, TY_CHAR) + call salloc (sf, SF, TY_STRUCT) + call salloc (sfd, SFD, TY_STRUCT) + call salloc (sfds, 1, TY_POINTER) + call aclri (Memi[sf], SF) + call aclri (Memi[sfd], SFD) + Memi[sfds] = sfd + + # Initialize parameters. + SF_TASK(sf) = PSFMEASURE + SF_WCODE(sf) = strdic ("FWHM", SF_WTYPE(sf), SF_SZWTYPE, SF_WTYPES) + SF_SCALE(sf) = 1. + SF_LEVEL(sf) = level + SF_BETA(sf) = beta + SF_RADIUS(sf) = radius + SF_SBUF(sf) = sbuffer + SF_SWIDTH(sf) = swidth + SF_SAT(sf) = saturation + SF_NIT(sf) = nit + SF_OVRPLT(sf) = NO + SF_NCOLS(sf) = IM_LEN(im,1) + SF_NLINES(sf) = IM_LEN(im,2) + SF_XF(sf) = (IM_LEN(im,1) + 1) / 2. + SF_YF(sf) = (IM_LEN(im,2) + 1) / 2. + ignore_sat = false + + call imstats (im, IM_IMAGENAME, SFD_IMAGE(sfd), SF_SZFNAME) + SFD_ID(sfd) = 1 + SFD_X(sfd) = xc + SFD_Y(sfd) = yc + SFD_F(sfd) = INDEF + SFD_STATUS(sfd) = 0 + SFD_SFS(sfd) = NULL + SFD_SFF(sfd) = NULL + SFD_SFI(sfd) = NULL + + if (SF_LEVEL(sf) > 1.) + SF_LEVEL(sf) = SF_LEVEL(sf) / 100. + SF_LEVEL(sf) = max (0.05, min (0.95, SF_LEVEL(sf))) + + # Evaluate PSF data. + iferr { + do i = 1, SF_NIT(sf) { + if (i == 1) + SFD_RADIUS(sfd) = SF_RADIUS(sf) + else + SFD_RADIUS(sfd) = 3. * SFD_DFWHM(sfd) + SFD_NPMAX(sfd) = stf_r2i (SFD_RADIUS(sfd)) + 1 + SFD_NP(sfd) = SFD_NPMAX(sfd) + call stf_find (sf, sfd, im) + call stf_bkgd (sf, sfd) + if (SFD_NSAT(sfd) > 0 && i == 1) { + if (ignore_sat) + call error (0, + "Saturated pixels found - ignoring object") + else + call eprintf ( + "WARNING: Saturated pixels found.\n") + } + call stf_profile (sf, sfd) + call stf_widths (sf, sfd) + call stf_fwhms (sf, sfd) + } + + # Set output results. + radius = SFD_RADIUS(sfd) + bkg = SFD_BKGD(sfd) + renclosed = SFD_R(sfd) + dfwhm = SFD_DFWHM(sfd) + mfwhm = SFD_MFWHM(sfd) + gfwhm = SFD_GFWHM(sfd) + + call asifree (SFD_ASI1(sfd)) + call asifree (SFD_ASI2(sfd)) + } then + call erract (EA_WARN) + + # Finish up + call stf_free (sf) + call sfree (sp) +end + + +# STF_FREE -- Free the starfocus data structures. + +procedure stf_free (sf) + +pointer sf #I Starfocus structure +int i + +begin + do i = 1, SF_NSTARS(sf) + call mfree (SF_SFS(sf,i), TY_STRUCT) + do i = 1, SF_NFOCUS(sf) + call mfree (SF_SFF(sf,i), TY_STRUCT) + do i = 1, SF_NIMAGES(sf) + call mfree (SF_SFI(sf,i), TY_STRUCT) + call mfree (SF_STARS(sf), TY_POINTER) + call mfree (SF_FOCUS(sf), TY_POINTER) + call mfree (SF_IMAGES(sf), TY_POINTER) + SF_NSTARS(sf) = 0 + SF_NFOCUS(sf) = 0 + SF_NIMAGES(sf) = 0 +end diff --git a/pkg/images/tv/imexamine/stfprofile.x b/pkg/images/tv/imexamine/stfprofile.x new file mode 100644 index 00000000..d26c085d --- /dev/null +++ b/pkg/images/tv/imexamine/stfprofile.x @@ -0,0 +1,1189 @@ +include <imhdr.h> +include <mach.h> +include <math.h> +include <math/iminterp.h> +include <math/nlfit.h> +include "starfocus.h" + + +# STF_FIND -- Find the object and return the data raster and object center. +# STF_BKGD -- Compute the background. +# STF_PROFILE -- Compute enclosed flux profile, derivative, and moments. +# STF_NORM -- Renormalized enclosed flux profile +# STF_WIDTHS -- Set widths. +# STF_I2R -- Radius from sample index. +# STF_R2I -- Sample index from radius. +# STF_R2N -- Number of subsamples from radius. +# STF_MODEL -- Return model values. +# STF_DFWHM -- Direct FWHM from profile. +# STF_FWHMS -- Measure FWHM vs level. +# STF_RADIUS -- Measure the radius at the specified level. +# STF_FIT -- Fit model. +# STF_GAUSS1 -- Gaussian function used in NLFIT. +# STF_GAUSS2 -- Gaussian function and derivatives used in NLFIT. +# STF_MOFFAT1 -- Moffat function used in NLFIT. +# STF_MOFFAT2 -- Moffat function and derivatives used in NLFIT. + + +# STF_FIND -- Find the object and return the data raster and object center. +# Centering uses centroid of marginal distributions of data above the mean. + +procedure stf_find (sf, sfd, im) + +pointer sf #I Starfocus pointer +pointer sfd #I Object pointer +pointer im #I Image pointer + +long lseed +int i, j, k, x1, x2, y1, y2, nx, ny, npts +real radius, buffer, width, xc, yc, xlast, ylast, r1, r2 +real mean, sum, sum1, sum2, sum3, asumr(), urand() +pointer data, ptr, imgs2r() +errchk imgs2r + +begin + radius = max (3., SFD_RADIUS(sfd)) + buffer = SF_SBUF(sf) + width = SF_SWIDTH(sf) + + xc = SFD_X(sfd) + yc = SFD_Y(sfd) + r1 = radius + buffer + width + r2 = radius + + # Iterate on the center finding. + do k = 1, 3 { + + # Extract region around current center. + xlast = xc + ylast = yc + + x1 = max (1-NBNDRYPIX, nint (xc - r2)) + x2 = min (IM_LEN(im,1)+NBNDRYPIX, nint (xc + r2)) + nx = x2 - x1 + 1 + y1 = max (1-NBNDRYPIX, nint (yc - r2)) + y2 = min (IM_LEN(im,2)+NBNDRYPIX, nint (yc + r2)) + ny = y2 - y1 + 1 + npts = nx * ny + data = imgs2r (im, x1, x2, y1, y2) + + # Find center of gravity of marginal distributions above mean. + npts = nx * ny + sum = asumr (Memr[data], npts) + mean = sum / nx + sum1 = 0. + sum2 = 0. + + do i = x1, x2 { + ptr = data + i - x1 + sum3 = 0. + do j = y1, y2 { + sum3 = sum3 + Memr[ptr] + ptr = ptr + nx + } + sum3 = sum3 - mean + if (sum3 > 0.) { + sum1 = sum1 + i * sum3 + sum2 = sum2 + sum3 + } + } + if (sum2 <= 0) + call error (1, "Centering failed to converge") + xc = sum1 / sum2 + if (xlast - xc > 0.2 * nx) + xc = xlast - 0.2 * nx + if (xc - xlast > 0.2 * nx) + xc = xlast + 0.2 * nx + + ptr = data + mean = sum / ny + sum1 = 0. + sum2 = 0. + do j = y1, y2 { + sum3 = 0. + do i = x1, x2 { + sum3 = sum3 + Memr[ptr] + ptr = ptr + 1 + } + sum3 = sum3 - mean + if (sum3 > 0.) { + sum1 = sum1 + j * sum3 + sum2 = sum2 + sum3 + } + } + if (sum2 <= 0) + call error (1, "Centering failed to converge") + yc = sum1 / sum2 + if (ylast - yc > 0.2 * ny) + yc = ylast - 0.2 * ny + if (yc - ylast > 0.2 * ny) + yc = ylast + 0.2 * ny + + if (nint(xc) == nint(xlast) && nint(yc) == nint(ylast)) + break + } + + # Get a new centered raster if necessary. + if (nint(xc) != nint(xlast) || nint(yc) != nint(ylast) || r2 < r1) { + x1 = max (1-NBNDRYPIX, nint (xc - r1)) + x2 = min (IM_LEN(im,1)+NBNDRYPIX, nint (xc + r1)) + nx = x2 - x1 + 1 + y1 = max (1-NBNDRYPIX, nint (yc - r1)) + y2 = min (IM_LEN(im,2)+NBNDRYPIX, nint (yc + r1)) + ny = y2 - y1 + 1 + npts = nx * ny + data = imgs2r (im, x1, x2, y1, y2) + } + + # Add a dither for integer data. The random numbers are always + # the same to provide reproducibility. + + i = IM_PIXTYPE(im) + if (i == TY_SHORT || i == TY_INT || i == TY_LONG) { + lseed = 1 + do i = 0, npts-1 + Memr[data+i] = Memr[data+i] + urand(lseed) - 0.5 + } + + SFD_DATA(sfd) = data + SFD_X1(sfd) = x1 + SFD_X2(sfd) = x2 + SFD_Y1(sfd) = y1 + SFD_Y2(sfd) = y2 + SFD_X(sfd) = xc + SFD_Y(sfd) = yc +end + + +# STF_BKGD -- Compute the background. +# A mode is estimated from the minimum slope in the sorted background pixels +# with a bin width of 5%. + +procedure stf_bkgd (sf, sfd) + +pointer sf #I Parameter structure +pointer sfd #I Star structure + +int i, j, x1, x2, y1, y2, xc, yc, nx, ny, npts, ns, nsat +real sat, bkgd, miso +real r, r1, r2, r3, dx, dy, dz +pointer sp, data, bdata, ptr + +begin + data = SFD_DATA(sfd) + x1 = SFD_X1(sfd) + x2 = SFD_X2(sfd) + y1 = SFD_Y1(sfd) + y2 = SFD_Y2(sfd) + xc = SFD_X(sfd) + yc = SFD_Y(sfd) + + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + npts = nx * ny + + ns = 0 + nsat = 0 + r1 = SFD_RADIUS(sfd) ** 2 + r2 = (SFD_RADIUS(sfd) + SF_SBUF(sf)) ** 2 + r3 = (SFD_RADIUS(sfd) + SF_SBUF(sf) + SF_SWIDTH(sf)) ** 2 + sat = SF_SAT(sf) + if (IS_INDEF(sat)) + sat = MAX_REAL + + call smark (sp) + call salloc (bdata, npts, TY_REAL) + + ptr = data + do j = y1, y2 { + dy = (yc - j) ** 2 + do i = x1, x2 { + dx = (xc - i) ** 2 + r = dx + dy + if (r <= r1) { + if (Memr[ptr] >= sat) + nsat = nsat + 1 + } else if (r >= r2 && r <= r3) { + Memr[bdata+ns] = Memr[ptr] + ns = ns + 1 + } + ptr = ptr + 1 + } + } + + if (ns > 9) { + call asrtr (Memr[bdata], Memr[bdata], ns) + r = Memr[bdata+ns-1] - Memr[bdata] + bkgd = Memr[bdata] + r / 2 + miso = r / 2 + + j = 1 + 0.50 * ns + do i = 0, ns - j { + dz = Memr[bdata+i+j-1] - Memr[bdata+i] + if (dz < r) { + r = dz + bkgd = Memr[bdata+i] + dz / 2 + miso = dz / 2 + } + } + } else { + bkgd = 0. + miso = 0. + } + + SFD_BKGD1(sfd) = bkgd + SFD_BKGD(sfd) = bkgd + SFD_MISO(sfd) = miso + SFD_NSAT(sfd) = nsat + + call sfree (sp) +end + + +# STF_PROFILE -- Compute enclosed flux profile, derivative, direct FWHM, and +# profile moments.. +# 1. The flux profile is normalized at the maximum value. +# 2. The radial profile is computed from the numerical derivative of the +# enclose flux profile. + +procedure stf_profile (sf, sfd) + +pointer sf #I Parameter structure +pointer sfd #I Star structure + +int np +real radius, xc, yc + +int i, j, k, l, m, ns, nx, ny, x1, x2, y1, y2 +real bkgd, miso, sigma, peak +real r, r1, r2, r3, dx, dy, dx1, dx2, dy1, dy2, dz, xx, yy, xy, ds, da +pointer sp, data, profile, ptr, asi, msi, gs +int stf_r2n() +real asieval(), msieval(), gseval(), stf_i2r(), stf_r2i() +errchk asiinit, asifit, msiinit, msifit, gsrestore + +real gsdata[24] +data gsdata/ 1., 4., 4., 1., 0., 0.6726812, 1., 2., 1.630641, 0.088787, + 0.00389378, -0.001457133, 0.3932125, -0.1267456, -0.004864541, + 0.00249941, 0.03078612, 0.02731274, -4.875850E-4, 2.307464E-4, + -0.002134843, 0.007603908, -0.002552385, -8.010564E-4/ + +begin + data = SFD_DATA(sfd) + x1 = SFD_X1(sfd) + x2 = SFD_X2(sfd) + y1 = SFD_Y1(sfd) + y2 = SFD_Y2(sfd) + xc = SFD_X(sfd) + yc = SFD_Y(sfd) + bkgd = SFD_BKGD(sfd) + miso = SFD_MISO(sfd) + radius = SFD_RADIUS(sfd) + np = SFD_NP(sfd) + + nx = x2 - x1 + 1 + ny = y2 - y1 + 1 + + # Use an image interpolator fit to the data. + call msiinit (msi, II_BISPLINE3) + call msifit (msi, Memr[data], nx, ny, nx) + + # To avoid trying to interpolate outside the center of the + # edge pixels, a requirement of the interpolator functions, + # we reset the data limits. + x1 = x1 + 1 + x2 = x2 - 1 + y1 = y1 + 1 + y2 = y2 - 1 + + # Compute the enclosed flux profile, its derivative, and moments. + call smark (sp) + call salloc (profile, np, TY_REAL) + call aclrr (Memr[profile], np) + + xx = 0. + yy = 0. + xy = 0. + do j = y1, y2 { + ptr = data + (j-y1+1)*nx + 1 + dy = j - yc + do i = x1, x2 { + dx = i - xc + + # Set the subpixel sampling which may be a function of radius. + r = sqrt (dx * dx + dy * dy) + ns = stf_r2n (r) + ds = 1. / ns + da = ds * ds + dz = 0.5 + 0.5 * ds + + # Sum the interpolator values over the subpixels and compute + # an offset to give the correct total for the pixel. + + r2 = 0. + dy1 = dy - dz + do l = 1, ns { + dy1 = dy1 + ds + dy2 = dy1 * dy1 + dx1 = dx - dz + do k = 1, ns { + dx1 = dx1 + ds + dx2 = dx1 * dx1 + r1 = msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2) + r2 = r2 + r1 + } + } + + r1 = Memr[ptr] - bkgd + ptr = ptr + 1 + r2 = r1 - r2 * da + + # Accumulate the enclosed flux over the sub pixels. + dy1 = dy - dz + do l = 1, ns { + dy1 = dy1 + ds + dy2 = dy1 * dy1 + dx1 = dx - dz + do k = 1, ns { + dx1 = dx1 + ds + dx2 = dx1 * dx1 + r = max (0., sqrt (dx2 + dy2) - ds / 2) + if (r < radius) { + r1 = da * (msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2) + + r2) + + # Use approximation for fractions of a subpixel. + for (m=stf_r2i(r)+1; m<=np; m=m+1) { + r3 = (stf_i2r (real(m)) - r) / ds + if (r3 >= 1.) + break + Memr[profile+m-1] = Memr[profile+m-1] + r3 * r1 + } + + # The subpixel is completely within these radii. + for (; m<=np; m=m+1) + Memr[profile+m-1] = Memr[profile+m-1] + r1 + + # Accumulate the moments above an isophote. + if (r1 > miso) { + xx = xx + dx2 * r1 + yy = yy + dy2 * r1 + xy = xy + dx1 * dy1 * r1 + } + } + } + } + } + } + + call msifree (msi) + + # Compute the ellipticity and position angle from the moments. + r = (xx + yy) + if (r > 0.) { + r1 = (xx - yy) / r + r2 = 2 * xy / r + SFD_E(sfd) = sqrt (r1**2 + r2**2) + SFD_PA(sfd) = RADTODEG (atan2 (r2, r1) / 2.) + } else { + SFD_E(sfd) = 0. + SFD_PA(sfd) = 0. + } + + # The magnitude and profile normalization is from the max enclosed flux. + call alimr (Memr[profile], np, r, SFD_M(sfd)) + if (SFD_M(sfd) <= 0.) + call error (1, "Invalid flux profile") + call adivkr (Memr[profile], SFD_M(sfd), Memr[profile], np) + + # Fit interpolator to the enclosed flux profile. + call asiinit (asi, II_SPLINE3) + call asifit (asi, Memr[profile], np) + SFD_ASI1(sfd) = asi + + # Estimate a gaussian sigma (actually sqrt(2)*sigma) and if it is + # it is small subtract the gaussian so that the image interpolator + # can more accurately estimate subpixel values. + + #call stf_radius (sf, sfd, SF_LEVEL(sf), r) + #sigma = r / sqrt (log (1/(1-SF_LEVEL(sf)))) + call stf_radius (sf, sfd, 0.8, r) + r = r / SF_SCALE(sf) + sigma = 2 * r * sqrt (log(2.) / log (1/(1-0.8))) + if (sigma < 5.) { + if (sigma <= 2.) { + call gsrestore (gs, gsdata) + dx = xc - nint (xc) + dy = yc - nint (yc) + r = sqrt (dx * dx + dy * dy) + dx = 1. + ds = abs (sigma - gseval (gs, r, dx)) + for (da = 1.; da <= 2.; da = da + .01) { + dz = abs (sigma - gseval (gs, r, da)) + if (dz < ds) { + ds = dz + dx = da + } + } + sigma = dx + call gsfree (gs) + } + + sigma = sigma / (2 * sqrt (log(2.))) + sigma = sigma * sigma + + # Compute the peak that gives the correct central pixel value. + i = nint (xc) + j = nint (yc) + dx = i - xc + dy = j - yc + r = sqrt (dx * dx + dy * dy) + ns = stf_r2n (r) + ds = 1. / ns + da = ds * ds + dz = 0.5 + 0.5 * ds + + r1 = 0. + dy1 = dy - dz + do l = 1, ns { + dy1 = dy1 + ds + dy2 = dy1 * dy1 + dx1 = dx - dz + do k = 1, ns { + dx1 = dx1 + ds + dx2 = dx1 * dx1 + r2 = (dx2 + dy2) / sigma + if (r2 < 25.) + r1 = r1 + exp (-r2) + } + } + ptr = data + (j - y1 + 1) * nx + (i - x1 + 1) + peak = (Memr[ptr] - bkgd) / (r1 * da) + + # Subtract the gaussian from the data. + do j = y1, y2 { + ptr = data + (j - y1 + 1) * nx + 1 + dy = j - yc + do i = x1, x2 { + dx = i - xc + r = sqrt (dx * dx + dy * dy) + ns = stf_r2n (r) + ds = 1. / ns + da = ds * ds + dz = 0.5 + 0.5 * ds + + r1 = 0. + dy1 = dy - dz + do l = 1, ns { + dy1 = dy1 + ds + dy2 = dy1 * dy1 + dx1 = dx - dz + do k = 1, ns { + dx1 = dx1 + ds + dx2 = dx1 * dx1 + r2 = (dx2 + dy2) / sigma + if (r2 < 25.) + r1 = r1 + peak * exp (-r2) + } + } + Memr[ptr] = Memr[ptr] - r1 * da + ptr = ptr + 1 + } + } + + # Fit the image interpolator to the residual data. + call msiinit (msi, II_BISPLINE3) + call msifit (msi, Memr[data], nx, ny, nx) + + # Recompute the enclosed flux profile and moments + # using the gaussian plus image interpolator fit to the residuals. + + call aclrr (Memr[profile], np) + + xx = 0. + yy = 0. + xy = 0. + do j = y1, y2 { + ptr = data + (j - y1 + 1) * nx + 1 + dy = j - yc + do i = x1, x2 { + dx = i - xc + r = sqrt (dx * dx + dy * dy) + ns = stf_r2n (r) + ds = 1. / ns + da = ds * ds + dz = 0.5 + 0.5 * ds + + # Compute interpolator correction. + r2 = 0. + dy1 = dy - dz + do l = 1, ns { + dy1 = dy1 + ds + dx1 = dx - dz + do k = 1, ns { + dx1 = dx1 + ds + r1 = msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2) + r2 = r2 + r1 + } + } + + r1 = Memr[ptr] - bkgd + ptr = ptr + 1 + r2 = r1 - r2 * da + + # Accumulate the enclosed flux and moments. + dy1 = dy - dz + do l = 1, ns { + dy1 = dy1 + ds + dy2 = dy1 * dy1 + dx1 = dx - dz + do k = 1, ns { + dx1 = dx1 + ds + dx2 = dx1 * dx1 + r3 = (dx2 + dy2) / sigma + if (r3 < 25.) + r3 = peak * exp (-r3) + else + r3 = 0. + r = max (0., sqrt (dx2 + dy2) - ds / 2) + if (r < radius) { + r1 = msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2) + r1 = da * (r1 + r2 + r3) + + for (m=stf_r2i(r)+1; m<=np; m=m+1) { + r3 = (stf_i2r (real(m)) - r) / ds + if (r3 >= 1.) + break + Memr[profile+m-1] = Memr[profile+m-1] + + r3 * r1 + } + for (; m<=np; m=m+1) + Memr[profile+m-1] = Memr[profile+m-1] + r1 + + if (r1 > miso) { + xx = xx + dx2 * r1 + yy = yy + dy2 * r1 + xy = xy + dx1 * dy1 * r1 + } + } + } + } + } + } + + call msifree (msi) + + # Recompute the moments, magnitude, normalized flux, and interp. + r = (xx + yy) + if (r > 0.) { + r1 = (xx - yy) / r + r2 = 2 * xy / r + SFD_E(sfd) = sqrt (r1**2 + r2**2) + SFD_PA(sfd) = RADTODEG (atan2 (r2, r1) / 2.) + } else { + SFD_E(sfd) = 0. + SFD_PA(sfd) = 0. + } + + call alimr (Memr[profile], np, r, SFD_M(sfd)) + if (SFD_M(sfd) <= 0.) + call error (1, "Invalid flux profile") + call adivkr (Memr[profile], SFD_M(sfd), Memr[profile], np) + + call asifit (asi, Memr[profile], np) + SFD_ASI1(sfd) = asi + } + + # Compute derivative of enclosed flux profile and fit an image + # interpolator. + + dx = 0.25 + Memr[profile] = 0. + ns = 0 + do i = 1, np { + r = stf_i2r (real(i)) + r2 = stf_r2i (r + dx) + if (r2 > np) { + k = i + break + } + r1 = stf_r2i (r - dx) + if (r1 < 1) { + if (i > 1) { + dy = asieval (asi, real(i)) / r**2 + Memr[profile] = (ns * Memr[profile] + dy) / (ns + 1) + ns = ns + 1 + } + j = i + } else { + dy = (asieval (asi, r2) - asieval (asi, r1)) / + (4 * r * dx) + Memr[profile+i-1] = dy + } + } + do i = 2, j + Memr[profile+i-1] = (Memr[profile+j] - Memr[profile]) / j * + (i - 1) + Memr[profile] + do i = k, np + Memr[profile+i-1] = Memr[profile+k-2] + + call adivkr (Memr[profile], SF_SCALE(sf)**2, Memr[profile], np) + call alimr (Memr[profile], np, SFD_YP1(sfd), SFD_YP2(sfd)) + call asiinit (asi, II_SPLINE3) + call asifit (asi, Memr[profile], np) + SFD_ASI2(sfd) = asi + #SF_XP1(sf) = j+1 + SF_XP1(sf) = 1 + SF_XP2(sf) = k-1 + + call sfree (sp) +end + + +# STF_NORM -- Renormalize the enclosed flux profile. + +procedure stf_norm (sf, sfd, x, y) + +pointer sf #I Parameter structure +pointer sfd #I Star structure +real x #I Radius +real y #I Flux + +int npmax, np +pointer asi + +int i, j, k +real r, r1, r2, dx, dy +pointer sp, profile +real asieval(), stf_i2r(), stf_r2i() +errchk asifit + +begin + npmax = SFD_NPMAX(sfd) + np = SFD_NP(sfd) + asi = SFD_ASI1(sfd) + + call smark (sp) + call salloc (profile, npmax, TY_REAL) + + # Renormalize the enclosed flux profile. + if (IS_INDEF(x) || x <= 0.) { + dy = SFD_BKGD(sfd) - SFD_BKGD1(sfd) + SFD_BKGD(sfd) = SFD_BKGD(sfd) - dy + do i = 1, npmax + Memr[profile+i-1] = asieval (asi, real(i)) + + dy * stf_i2r(real(i)) ** 2 + call alimr (Memr[profile], np, r1, r2) + call adivkr (Memr[profile], r2, Memr[profile], npmax) + } else if (IS_INDEF(y)) { + r = max (1., min (real(np), stf_r2i (x))) + r2 = asieval (asi, r) + if (r2 <= 0.) + return + do i = 1, npmax + Memr[profile+i-1] = asieval (asi, real(i)) + call adivkr (Memr[profile], r2, Memr[profile], npmax) + } else { + r = max (1., min (real(np), stf_r2i (x))) + r1 = asieval (asi, r) + dy = (y - r1) / x ** 2 + SFD_BKGD(sfd) = SFD_BKGD(sfd) - dy + do i = 1, npmax + Memr[profile+i-1] = asieval (asi, real(i)) + + dy * stf_i2r(real(i)) ** 2 + } + + call asifit (asi, Memr[profile], npmax) + SFD_ASI1(sfd) = asi + + # Compute derivative of enclosed flux profile and fit an image + # interpolator. + + dx = 0.25 + do i = 1, npmax { + r = stf_i2r (real(i)) + r2 = stf_r2i (r + dx) + if (r2 > np) { + k = i + break + } + r1 = stf_r2i (r - dx) + if (r1 < 1) { + if (i > 1) { + dy = asieval (asi, real(i)) / r**2 + Memr[profile] = dy + } + j = i + } else { + dy = (asieval (asi, r2) - asieval (asi, r1)) / + (4 * r * dx) + Memr[profile+i-1] = dy + } + } + do i = 2, j + Memr[profile+i-1] = (Memr[profile+j] - Memr[profile]) / j * + (i - 1) + Memr[profile] + do i = k, npmax + Memr[profile+i-1] = Memr[profile+k-2] + + call adivkr (Memr[profile], SF_SCALE(sf)**2, Memr[profile], np) + call alimr (Memr[profile], np, SFD_YP1(sfd), SFD_YP2(sfd)) + asi = SFD_ASI2(sfd) + call asifit (asi, Memr[profile], np) + SFD_ASI2(sfd) = asi + #SF_XP1(sf) = min (j+1, np) + SF_XP1(sf) = 1 + SF_XP2(sf) = min (k-1, np) + + call sfree (sp) +end + + +# STF_WIDTHS -- Set the widhts. + +procedure stf_widths (sf, sfd) + +pointer sf #I Main data structure +pointer sfd #I Star data structure + +errchk stf_radius, stf_dfwhm, stf_fit + +begin + call stf_radius (sf, sfd, SF_LEVEL(sf), SFD_R(sfd)) + call stf_dfwhm (sf, sfd) + call stf_fit (sf, sfd) + + switch (SF_WCODE(sf)) { + case 1: + SFD_W(sfd) = SFD_R(sfd) + case 2: + SFD_W(sfd) = SFD_DFWHM(sfd) + case 3: + SFD_W(sfd) = SFD_GFWHM(sfd) + case 4: + SFD_W(sfd) = SFD_MFWHM(sfd) + } +end + + +# STF_I2R -- Compute radius from sample index. + +real procedure stf_i2r (i) + +real i #I Index +real r #O Radius + +begin + if (i < 20) + r = 0.05 * i + else if (i < 30) + r = 0.1 * i - 1 + else if (i < 40) + r = 0.2 * i - 4 + else if (i < 50) + r = 0.5 * i - 16 + else + r = i - 41 + return (r) +end + + +# STF_R2I -- Compute sample index from radius. + +real procedure stf_r2i (r) + +real r #I Radius +real i #O Index + +begin + if (r < 1) + i = 20 * r + else if (r < 2) + i = 10 * (r + 1) + else if (r < 4) + i = 5 * (r + 4) + else if (r < 9) + i = 2 * (r + 16) + else + i = r + 41 + return (i) +end + + +# STF_R2N -- Compute number of subsamples from radius. + +int procedure stf_r2n (r) + +real r #I Radius +int n #O Number of subsamples + +begin + if (r < 1) + n = 20 + else if (r < 2) + n = 10 + else if (r < 4) + n = 5 + else if (r < 9) + n = 2 + else + n = 1 + return (n) +end + + +# STF_MODEL -- Return model value. + +procedure stf_model (sf, sfd, r, profile, flux) + +pointer sf #I Main data structure +pointer sfd #I Star data structure +real r #I Radius at level +real profile #I Profile value +real flux #I Enclosed flux value + +real x, x1, x2, r1, r2, dr + +begin + dr = 0.25 * SF_SCALE(sf) + r1 = r - dr + r2 = r + dr + if (r1 < 0.) { + r1 = dr + r2 = r1 + dr + } + + switch (SF_WCODE(sf)) { + case 3: + x = r**2 / (2. * SFD_SIGMA(sfd)**2) + if (x < 20.) + flux = 1 - exp (-x) + else + flux = 0. + + x1 = r1**2 / (2. * SFD_SIGMA(sfd)**2) + x2 = r2**2 / (2. * SFD_SIGMA(sfd)**2) + if (x2 < 20.) { + x1 = 1 - exp (-x1) + x2 = 1 - exp (-x2) + } else { + x1 = 1. + x2 = 1. + } + if (r <= dr) { + x1 = x1 / dr ** 2 + x2 = x2 / (4 * dr ** 2) + profile = (x2 - x1) / dr * r + x1 + } else { + profile = (x2 - x1) / (4 * r * dr) + } + default: + x = 1 + (r / SFD_ALPHA(sfd)) ** 2 + flux = 1 - x ** (1 - SFD_BETA(sfd)) + + x1 = 1 + (r1 / SFD_ALPHA(sfd)) ** 2 + x2 = 1 + (r2 / SFD_ALPHA(sfd)) ** 2 + x1 = 1 - x1 ** (1 - SFD_BETA(sfd)) + x2 = 1 - x2 ** (1 - SFD_BETA(sfd)) + if (r <= dr) { + x1 = x1 / dr ** 2 + x2 = x2 / (4 * dr ** 2) + profile = (x2 - x1) / dr * r + x1 + } else { + profile = (x2 - x1) / (4 * r * dr) + } + } +end + + +# STF_DFWHM -- Direct FWHM from profile. + +procedure stf_dfwhm (sf, sfd) + +pointer sf #I Main data structure +pointer sfd #I Star data structure + +int np +real r, rpeak, profile, peak, asieval(), stf_i2r() +pointer asi + +begin + asi = SFD_ASI2(sfd) + np = SFD_NP(sfd) + + rpeak = 1. + peak = 0. + for (r=1.; r <= np; r = r + 0.01) { + profile = asieval (asi, r) + if (profile > peak) { + rpeak = r + peak = profile + } + } + + peak = peak / 2. + for (r=rpeak; r <= np && asieval (asi, r) > peak; r = r + 0.01) + ; + + SFD_DFWHM(sfd) = 2 * stf_i2r (r) * SF_SCALE(sf) +end + + +# STF_FWHMS -- Measure FWHM vs level. + +procedure stf_fwhms (sf, sfd) + +pointer sf #I Main data structure +pointer sfd #I Star data structure + +int i +real level, r + +begin + do i = 1, 19 { + level = i * 0.05 + call stf_radius (sf, sfd, level, r) + switch (SF_WCODE(sf)) { + case 3: + SFD_FWHM(sfd,i) = 2 * r * sqrt (log (2.) / log (1/(1-level))) + default: + r = r / sqrt ((1.-level)**(1./(1.-SFD_BETA(sfd))) - 1.) + SFD_FWHM(sfd,i) = 2 * r * sqrt (2.**(1./SFD_BETA(sfd))-1.) + } + } +end + + +# STF_RADIUS -- Measure the radius at the specified level. + +procedure stf_radius (sf, sfd, level, r) + +pointer sf #I Main data structure +pointer sfd #I Star data structure +real level #I Level to measure +real r #O Radius + +int np +pointer asi +real f, fmax, rmax, asieval(), stf_i2r() + +begin + np = SFD_NP(sfd) + asi = SFD_ASI1(sfd) + + for (r=1; r <= np && asieval (asi, r) < level; r = r + 0.01) + ; + if (r > np) { + fmax = 0. + rmax = 0. + for (r=1; r <= np; r = r + 0.01) { + f = asieval (asi, r) + if (f > fmax) { + fmax = f + rmax = r + } + } + r = rmax + } + r = stf_i2r (r) * SF_SCALE(sf) +end + + +# STF_FIT -- Fit models to enclosed flux. + +procedure stf_fit (sf, sfd) + +pointer sf #I Main data structure +pointer sfd #I Star data structure + +int i, j, n, np, pfit[2] +real beta, z, params[3] +pointer asi, nl +pointer sp, x, y, w + +int locpr() +real asieval(), stf_i2r() +extern stf_gauss1(), stf_gauss2(), stf_moffat1(), stf_moffat2() +errchk nlinitr, nlfitr + +data pfit/2,3/ + +begin + np = SFD_NP(sfd) + asi = SFD_ASI1(sfd) + + call smark (sp) + call salloc (x, np, TY_REAL) + call salloc (y, np, TY_REAL) + call salloc (w, np, TY_REAL) + + n = 0 + j = 0 + do i = 1, np { + z = 1. - max (0., asieval (asi, real(i))) + if (n > np/3 && z < 0.5) + break + if ((n < np/3 && z > 0.01) || z > 0.5) + j = n + + Memr[x+n] = stf_i2r (real(i)) * SF_SCALE(sf) + Memr[y+n] = z + Memr[w+n] = 1. + n = n + 1 + } + + # Gaussian. + np = 1 + params[2] = Memr[x+j] / sqrt (2. * log (1./min(0.99,Memr[y+j]))) + params[1] = 1 + call nlinitr (nl, locpr (stf_gauss1), locpr (stf_gauss2), + params, params, 2, pfit, np, .001, 100) + call nlfitr (nl, Memr[x], Memr[y], Memr[w], n, 1, WTS_USER, i) + if (i != SINGULAR && i != NO_DEG_FREEDOM) { + call nlpgetr (nl, params, i) + if (params[2] < 0.) + params[2] = Memr[x+j] / sqrt (2. * log (1./min(0.99,Memr[y+j]))) + } + SFD_SIGMA(sfd) = params[2] + SFD_GFWHM(sfd) = 2 * SFD_SIGMA(sfd) * sqrt (2. * log (2.)) + + # Moffat. + if (SF_BETA(sf) < 1.1) { + call nlfreer (nl) + call sfree (sp) + call error (1, "Cannot measure FWHM - Moffat beta too small") + } + + beta = SF_BETA(sf) + if (IS_INDEFR(beta)) { + beta = 2.5 + np = 2 + } else { + np = 1 + } + params[3] = 1 - beta + params[2] = Memr[x+j] / sqrt (min(0.99,Memr[y+j])**(1./params[3]) - 1.) + params[1] = 1 + call nlinitr (nl, locpr (stf_moffat1), locpr (stf_moffat2), + params, params, 3, pfit, np, .001, 100) + call nlfitr (nl, Memr[x], Memr[y], Memr[w], n, 1, WTS_USER, i) + if (i != SINGULAR && i != NO_DEG_FREEDOM) { + call nlpgetr (nl, params, i) + if (params[2] < 0.) { + params[3] = 1. - beta + params[2] = Memr[x+j] / + sqrt (min(0.99,Memr[y+j])**(1./params[3]) - 1.) + } + } + SFD_ALPHA(sfd) = params[2] + SFD_BETA(sfd) = 1 - params[3] + SFD_MFWHM(sfd) = 2 * SFD_ALPHA(sfd) * sqrt (2.**(1./SFD_BETA(sfd))-1.) + + call nlfreer (nl) + call sfree (sp) +end + + +# STF_GAUSS1 -- Gaussian function used in NLFIT. The parameters are the +# amplitude and sigma and the input variable is the radius. + +procedure stf_gauss1 (x, nvars, p, np, z) + +real x[nvars] #I Input variables +int nvars #I Number of variables +real p[np] #I Parameter vector +int np #I Number of parameters +real z #O Function return + +real r2 + +begin + r2 = x[1]**2 / (2 * p[2]**2) + if (abs (r2) > 20.) + z = 0. + else + z = p[1] * exp (-r2) +end + + +# STF_GAUSS2 -- Gaussian function and derivatives used in NLFIT. The parameters +# are the amplitude and sigma and the input variable is the radius. + +procedure stf_gauss2 (x, nvars, p, dp, np, z, der) + +real x[nvars] #I Input variables +int nvars #I Number of variables +real p[np] #I Parameter vector +real dp[np] #I Dummy array of parameters increments +int np #I Number of parameters +real z #O Function return +real der[np] #O Derivatives + +real r2 + +begin + r2 = x[1]**2 / (2 * p[2]**2) + if (abs (r2) > 20.) { + z = 0. + der[1] = 0. + der[2] = 0. + } else { + der[1] = exp (-r2) + z = p[1] * der[1] + der[2] = z * 2 * r2 / p[2] + } +end + + +# STF_MOFFAT1 -- Moffat function used in NLFIT. The parameters are the +# amplitude, alpha squared, and beta and the input variable is the radius. + +procedure stf_moffat1 (x, nvars, p, np, z) + +real x[nvars] #I Input variables +int nvars #I Number of variables +real p[np] #I Parameter vector +int np #I Number of parameters +real z #O Function return + +real y + +begin + y = 1 + (x[1] / p[2]) ** 2 + if (abs (y) > 20.) + z = 0. + else + z = p[1] * y ** p[3] +end + + +# STF_MOFFAT2 -- Moffat function and derivatives used in NLFIT. The +# parameters are the amplitude, alpha squared, and beta and the input +# variable is the radius. + +procedure stf_moffat2 (x, nvars, p, dp, np, z, der) + +real x[nvars] #I Input variables +int nvars #I Number of variables +real p[np] #I Parameter vector +real dp[np] #I Dummy array of parameters increments +int np #I Number of parameters +real z #O Function return +real der[np] #O Derivatives + +real y + +begin + y = 1 + (x[1] / p[2]) ** 2 + if (abs (y) > 20.) { + z = 0. + der[1] = 0. + der[2] = 0. + der[3] = 0. + } else { + der[1] = y ** p[3] + z = p[1] * der[1] + der[2] = -2 * z / y * p[3] / p[2] * (x[1] / p[2]) ** 2 + der[3] = z * log (y) + } +end diff --git a/pkg/images/tv/imexamine/t_imexam.x b/pkg/images/tv/imexamine/t_imexam.x new file mode 100644 index 00000000..089e74fc --- /dev/null +++ b/pkg/images/tv/imexamine/t_imexam.x @@ -0,0 +1,352 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <gset.h> +include <imhdr.h> +include "imexam.h" + +define HELP "iraf$lib/scr/imexamine.key" +define PROMPT "imexamine options" +define SZ_IMLIST 512 + + +# T_IMEXAMINE -- Examine images using image display, graphics, and text output. + +procedure t_imexamine () + +real x, y +pointer sp, cmd, imname, imlist, gp, ie, im +int curtype, key, redraw, mode, nframes, nargs + +bool clgetb() +pointer gopen(), ie_gimage(), imtopen() +int imd_wcsver(), ie_gcur(), ie_getnframes() +int btoi(), clgeti(), imtlen() + +begin + call smark (sp) + call salloc (ie, IE_LEN, TY_STRUCT) + call salloc (cmd, SZ_LINE, TY_CHAR) + call salloc (imname, SZ_FNAME, TY_CHAR) + call salloc (imlist, SZ_IMLIST, TY_CHAR) + + # Initialize the imexamine descriptor. + call aclri (Memi[ie], IE_LEN) + + # Determine if we will be accessing the image display, and if so, + # the maximum number of frames to be accessed. + + IE_USEDISPLAY(ie) = btoi (clgetb ("use_display")) + if (IE_USEDISPLAY(ie) == YES) { + if (imd_wcsver() == 0) + ; + iferr (nframes = ie_getnframes (ie)) { + call eprintf ("cannot access display\n") + IE_USEDISPLAY(ie) = NO + } + } + + # Get the list of images to be examined, if given on the command + # line. If no images are explicitly listed use the display to + # determine the images to be examined. + + nargs = clgeti ("$nargs") + if (nargs > 0) { + call clgstr ("input", Memc[imlist], SZ_IMLIST) + IE_LIST(ie) = imtopen (Memc[imlist]) + IE_LISTLEN(ie) = imtlen (IE_LIST(ie)) + IE_INDEX(ie) = 1 + + if (nargs >= 1) { + # Set user specified display frame. + IE_DFRAME(ie) = 100 * clgeti ("frame") + 1 + IE_NEWFRAME(ie) = IE_DFRAME(ie) + if (IE_USEDISPLAY(ie) == YES) { + nframes = max (IE_NEWFRAME(ie)/100, nframes) + IE_NFRAMES(ie) = nframes + } + } else { + # If we have to display an image and no frame was specified, + # default to frame 1 (should use the current display frame + # but we don't have a cursor read yet to tell us what it is). + + IE_DFRAME(ie) = 101 + IE_NEWFRAME(ie) = 101 + } + + } else { + IE_INDEX(ie) = 1 + IE_DFRAME(ie) = 101 + IE_NEWFRAME(ie) = 101 + } + + # Set the wcs, logfile and graphics. + call clgstr ("wcs", IE_WCSNAME(ie), IE_SZFNAME) + IE_LOGFD(ie) = NULL + call clgstr ("logfile", IE_LOGFILE(ie), IE_SZFNAME) + if (clgetb ("keeplog")) + iferr (call ie_openlog (ie)) + call erract (EA_WARN) + + call clgstr ("graphics", Memc[cmd], SZ_LINE) + gp = gopen (Memc[cmd], NEW_FILE+AW_DEFER, STDGRAPH) + + # Initialize the data structure. + IE_IM(ie) = NULL + IE_DS(ie) = NULL + IE_PP(ie) = NULL + IE_MAPFRAME(ie) = 0 + IE_NFRAMES(ie) = nframes + IE_ALLFRAMES(ie) = btoi (clgetb ("allframes")) + IE_GTYPE(ie) = NULL + IE_XORIGIN(ie) = 0. + IE_YORIGIN(ie) = 0. + + # Access the first image. If an image list was specified and the + # display is being used, this will set the display frame to the first + # image listed, or display the first image if not already loaded into + # the display. + + if (IE_LIST(ie) != NULL) + im = ie_gimage (ie, YES) + + # Enter the cursor loop. The commands are returned by the + # IE_GCUR procedure. + + x = 1. + y = 1. + redraw = NO + curtype = 'i' + mode = NEW_FILE + + while (ie_gcur (ie, curtype, x,y, key, Memc[cmd], SZ_LINE) != EOF) { + # Check to see if the user has changed frames on us while in + # examine-image-list mode. + + if (IE_USEDISPLAY(ie) == YES && IE_LIST(ie) != NULL && + IE_NEWFRAME(ie)/100 != IE_MAPFRAME(ie)/100) { + call ie_imname (IE_DS(ie), IE_NEWFRAME(ie), Memc[imname], + SZ_FNAME) + call ie_addimage (ie, Memc[imname], imlist) + } + + # Set workstation state. + switch (key) { + case 'a', 'b', 'd', 'm', 't', 'w', 'x', 'y', 'z', ',': + call gdeactivate (gp, 0) + } + + # Act on the command key. + switch (key) { + case '?': # Print help + call gpagefile (gp, HELP, PROMPT) + case ':': # Process colon commands + call ie_colon (ie, Memc[cmd], gp, redraw) + if (redraw == YES) { + x = INDEF + y = INDEF + } + case 'f': # Redraw frame + redraw = YES + x = INDEF + y = INDEF + case 'a': # Aperture photometry + call ie_rimexam (NULL, NULL, ie, x, y) + case ',': # Aperture photometry + call ie_qrimexam (NULL, NULL, ie, x, y) + + case 'b': # Print image region coordinates + call printf ("%4d %4d %4d %4d\n") + call pargi (IE_IX1(ie)) + call pargi (IE_IX2(ie)) + call pargi (IE_IY1(ie)) + call pargi (IE_IY2(ie)) + + if (IE_LOGFD(ie) != NULL) { + call fprintf (IE_LOGFD(ie), "%4d %4d %4d %4d\n") + call pargi (IE_IX1(ie)) + call pargi (IE_IX2(ie)) + call pargi (IE_IY1(ie)) + call pargi (IE_IY2(ie)) + } + + case 'c','e','h','j','k','s','l','r','u','v','.': # Graphs + IE_GTYPE(ie) = key + redraw = YES + + case 'd': # Load the display. + # Query the user for the frame to be loaded, the current + # display frame being the default. + + call clgstr ("image", Memc[imname], SZ_FNAME) + call clputi ("frame", IE_NEWFRAME(ie)/100) + IE_DFRAME(ie) = 100 * clgeti ("frame") + 1 + IE_NEWFRAME(ie) = IE_DFRAME(ie) + + if (IE_LIST(ie) != NULL) + call ie_addimage (ie, Memc[imname], imlist) + else + call ie_display (ie, Memc[imname], IE_DFRAME(ie)/100) + + case 'g': # Graphics cursor + curtype = 'g' + case 'i': # Image cursor + curtype = 'i' + case 'm': # Image statistics + call ie_statistics (ie, x, y) + + case 'n': # Next frame + if (IE_LIST(ie) != NULL) { + IE_INDEX(ie) = IE_INDEX(ie) + 1 + if (IE_INDEX(ie) > IE_LISTLEN(ie)) + IE_INDEX(ie) = 1 + } else { + IE_NEWFRAME(ie) = 100 * (IE_NEWFRAME(ie)/100 + 1) + 1 + if (IE_NEWFRAME(ie)/100 > IE_NFRAMES(ie)) + IE_NEWFRAME(ie) = 101 + } + im = ie_gimage (ie, YES) + + case 'o': # Overplot + mode = APPEND + + case 'p': # Previous frame + if (IE_LIST(ie) != NULL) { + IE_INDEX(ie) = IE_INDEX(ie) - 1 + if (IE_INDEX(ie) <= 0) + IE_INDEX(ie) = IE_LISTLEN(ie) + } else { + IE_NEWFRAME(ie) = 100 * (IE_NEWFRAME(ie)/100 - 1) + 1 + if (IE_NEWFRAME(ie)/100 <= 0) + IE_NEWFRAME(ie) = 100 * IE_NFRAMES(ie) + 1 + } + im = ie_gimage (ie, YES) + + case 'q': # Quit + break + + case 't': # Extract a section. + call ie_timexam (ie, x, y) + + case 'w': # Toggle logfile + if (IE_LOGFD(ie) == NULL) { + if (IE_LOGFILE(ie) == EOS) + call printf ("No log file defined\n") + else { + iferr (call ie_openlog (ie)) + call erract (EA_WARN) + } + } else { + call close (IE_LOGFD(ie)) + IE_LOGFD(ie) = NULL + call printf ("Logfile %s closed\n") + call pargstr (IE_LOGFILE(ie)) + } + + case 'x', 'y': # Positions + call ie_pos (ie, x, y, key) + case 'z': # Print grid + call ie_print (ie, x, y) + case 'I': # Immediate interrupt + call fatal (1, "Interrupt") + default: # Unrecognized command + call printf ("\007") + } + + switch (key) { + case '?', 'a', 'b', 'd', 'm', 'w', 'x', 'y', 'z', ',': + IE_LASTKEY(ie) = key + } + + # Draw or overplot a graph. + if (redraw == YES) { + switch (IE_GTYPE(ie)) { + case 'c': # column plot + call ie_cimexam (gp, mode, ie, x) + case 'e': # contour plot + call ie_eimexam (gp, mode, ie, x, y) + case 'h': # histogram plot + call ie_himexam (gp, mode, ie, x, y) + case 'j': # line plot + call ie_jimexam (gp, mode, ie, x, y, 1) + case 'k': # line plot + call ie_jimexam (gp, mode, ie, x, y, 2) + case 'l': # line plot + call ie_limexam (gp, mode, ie, y) + case 'r': # radial profile plot + call ie_rimexam (gp, mode, ie, x, y) + case 's': # surface plot + call ie_simexam (gp, mode, ie, x, y) + case 'u', 'v': # vector cut plot + call ie_vimexam (gp, mode, ie, x, y, IE_GTYPE(ie)) + case '.': # radial profile plot + call ie_qrimexam (gp, mode, ie, x, y) + } + redraw = NO + mode = NEW_FILE + } + } + + # Finish up. + call gclose (gp) + if (IE_IM(ie) != NULL && IE_IM(ie) != IE_DS(ie)) + call imunmap (IE_IM(ie)) + if (IE_MW(ie) != NULL) + call mw_close (IE_MW(ie)) + if (IE_PP(ie) != NULL) + call clcpset (IE_PP(ie)) + if (IE_DS(ie) != NULL) + call imunmap (IE_DS(ie)) + if (IE_LOGFD(ie) != NULL) + call close (IE_LOGFD(ie)) + if (IE_LIST(ie) != NULL) + call imtclose (IE_LIST(ie)) + call sfree (sp) +end + + +# IE_ADDIMAGE -- Add an image to the image list if not already present in +# the list, and display the image. + +procedure ie_addimage (ie, image, imlist) + +pointer ie #I imexamine descriptor +char image[ARB] #I image name +pointer imlist #I image list + +int i +bool inlist +pointer im, sp, lname +pointer ie_gimage(), imtopen() +int imtrgetim(), imtlen() +bool streq() + +begin + call smark (sp) + call salloc (lname, SZ_FNAME, TY_CHAR) + + # Is image already in list? + inlist = false + do i = 1, IE_LISTLEN(ie) { + if (imtrgetim (IE_LIST(ie), i, Memc[lname], SZ_FNAME) > 0) + if (streq (Memc[lname], image)) { + inlist = true + IE_INDEX(ie) = i + break + } + } + + # Add to list if missing. + if (!inlist) { + call strcat (",", Memc[imlist], SZ_IMLIST) + call strcat (image, Memc[imlist], SZ_IMLIST) + call imtclose (IE_LIST(ie)) + IE_LIST(ie) = imtopen (Memc[imlist]) + IE_LISTLEN(ie) = imtlen (IE_LIST(ie)) + IE_INDEX(ie) = IE_LISTLEN(ie) + } + + # Display the image. + im = ie_gimage (ie, YES) + call sfree (sp) +end diff --git a/pkg/images/tv/imexamine/x_imexam.x b/pkg/images/tv/imexamine/x_imexam.x new file mode 100644 index 00000000..100a6756 --- /dev/null +++ b/pkg/images/tv/imexamine/x_imexam.x @@ -0,0 +1 @@ +task imexamine = t_imexamine |