diff options
Diffstat (limited to 'noao/rv/rvidlines')
37 files changed, 7783 insertions, 0 deletions
diff --git a/noao/rv/rvidlines/Revisions b/noao/rv/rvidlines/Revisions new file mode 100644 index 00000000..9f93511a --- /dev/null +++ b/noao/rv/rvidlines/Revisions @@ -0,0 +1,52 @@ +idinit.x +reidentify.x + Changed aclrr to aclri and updated the 'b' key in reidentify.x. + (12/1/93, Valdes) + +idlog.x + Added the spectrum title and heliocentric Julian date to the log output. + (11/17/93, Valdes) + +iddblend.x + + Needed copy of ONEDSPEC deblend with names changes to protect the + innocent and avoid conflict with the version in RV. (8/25/93, Valdes) + +t_identify.x +t_reidentify.x +idcenter.x +idcolon.x +iddb.x +iddofit.x +iddoshift.x +idgdata.x +idgraph.x +ididentify.x +idinit.x +idlinelist.x +idlog.x +idmark.x +idrms.x +idshift.x +idshow.x +idvelocity.x + +idvhelio.x +peaks.gx +reidentify.x +mkpkg +identify.h +identify.par +reidentify.par +rvidlines.par + +rvreidlines.par + +identify.key +rvidlines.key + +doc/rvidlines.hlp + +doc/rvreidlines.hlp + + 1. A new cursor key was added for marking and deblending features. + 2. Two new feature types (ftype) were added: + gabsorption Absorption features fit by gaussian + gemission Emission features fit by gaussian + 3. The 'u' key was modified to allow setting a feature weight + 4. Two new tasks were added RVIDLINES and RVREIDLINES. These + are implemented as special entry points and a task flag. + (8/23/93, Valdes) diff --git a/noao/rv/rvidlines/idcenter.x b/noao/rv/rvidlines/idcenter.x new file mode 100644 index 00000000..05c636cb --- /dev/null +++ b/noao/rv/rvidlines/idcenter.x @@ -0,0 +1,257 @@ +include <gset.h> +include <smw.h> +include "identify.h" + +# ID_CENTER -- Locate the center of a feature. + +double procedure id_center (id, x, n, width, type, interactive) + +pointer id # ID pointer +double x[n] # Initial guess +int n # Number of features +real width # Feature width +int type # Feature type +int interactive # Interactive? + +int np1 +real value + +real center1d() +double smw_c1trand() + +begin + switch (type) { + case EMISSION, ABSORPTION: + np1 = NP1(ID_SH(id)) - 1 + value = smw_c1trand (ID_PL(id), x[1]) - np1 + value = center1d (value, IMDATA(id,1), ID_NPTS(id), + width, type, ID_CRADIUS(id), ID_THRESHOLD(id)) + if (IS_INDEF(value)) + return (INDEFD) + else + return (smw_c1trand (ID_LP(id), double(value+np1))) + case GEMISSION, GABSORPTION: + iferr (call id_gcenter (id, x, n, INDEF, INDEF, INDEF, INDEF, + width, type, interactive)) + return (INDEFD) + return (x[1]) + } +end + + +# ID_GCENTER -- Locate the center of a feature. + +procedure id_gcenter (id, x, ng, xa, ya, xb, yb, width, type, interactive) + +pointer id # ID pointer +double x[ng] # Initial guess +int ng # Number of features +real xa, ya, xb, yb # Background points +real width # Feature width +int type # Feature type +int interactive # Draw gaussian fit? + +int i, np1, x1, x2 +real ag, bg, pix, w, y +pointer sp, xr, xg, yg, sg, gp + +bool fp_equalr() +real id_model() +double smw_c1trand(), id_fitpt() + +errchk gcenter1d + +begin + call smark (sp) + call salloc (xr, ng, TY_REAL) + call salloc (xg, ng, TY_REAL) + call salloc (yg, ng, TY_REAL) + call salloc (sg, ng, TY_REAL) + + np1 = NP1(ID_SH(id)) - 1 + + # Compute background in logical units. + if (IS_INDEF(xa) || IS_INDEF(ya) || IS_INDEF(xb) || IS_INDEF(yb)) { + ag = INDEF + bg = INDEF + } else { + ag = smw_c1trand (ID_PL(id), double(xa)) - np1 + bg = smw_c1trand (ID_PL(id), double(xb)) - np1 + if (!fp_equalr (ag, bg)) { + bg = (yb - ya) / (bg - ag) + ag = ya - bg * ag + } else { + ag = INDEF + bg = INDEF + } + } + + do i = 1, ng + Memr[xr+i-1] = smw_c1trand (ID_PL(id), x[i]) - np1 + call gcenter1d (Memr[xr], ng, IMDATA(id,1), ID_NPTS(id), + width, type, ID_CRADIUS(id), ID_THRESHOLD(id), + x1, x2, Memr[xg], Memr[yg], Memr[sg], ag, bg) + do i = 1, ng + x[i] = smw_c1trand (ID_LP(id), double(Memr[xg+i-1]+np1)) + + if (interactive == YES) { + gp = ID_GP(id) + call gseti (gp, G_PLTYPE, 2) + call gseti (gp, G_PLCOLOR, 2) + pix = x1 + w = id_fitpt (id, smw_c1trand (ID_LP(id), double (pix+np1))) + y = id_model (pix, Memr[xg], Memr[yg], Memr[sg], ng) + ag + + bg * pix + call gamove (gp, w, y) + for (pix = x1; pix <= x2; pix = pix + .1) { + w = id_fitpt (id, smw_c1trand (ID_LP(id), double (pix+np1))) + y = id_model (pix, Memr[xg], Memr[yg], Memr[sg], ng) + + ag + bg * pix + call gadraw (gp, w, y) + } + call gseti (gp, G_PLTYPE, 3) + call gseti (gp, G_PLCOLOR, 3) + pix = x1 + w = id_fitpt (id, smw_c1trand (ID_LP(id), double (pix+np1))) + y = ag + bg * pix + call gamove (gp, w, y) + pix = x2 + w = id_fitpt (id, smw_c1trand (ID_LP(id), double (pix+np1))) + y = ag + bg * pix + call gadraw (gp, w, y) + call gseti (gp, G_PLTYPE, 1) + call gseti (gp, G_PLCOLOR, 1) + call gflush (gp) + } + + call sfree (sp) +end + + +define MIN_WIDTH 3. # Minimum centering width + + +# GCENTER1D -- Locate the center of a one dimensional feature by guassian fit. +# A value of INDEF is returned in the centering fails for any reason. +# This procedure just sets up the data and adjusts for emission or +# absorption features. The actual centering is done by GFIT. +# If width <= 1 return the nearest minima or maxima. + +procedure gcenter1d (x, ng, data, npts, width, type, radius, threshold, + x1, x2, xg, yg, sg, ag, bg) + +real x[ng] # Initial guess +int ng # Number of gaussians +real data[npts] # Data points +int npts # Number of data points +real width # Feature width +int type # Feature type +real radius # Centering radius +real threshold # Minimum range in feature +int x1, x2 # Fitting region +real xg[ng], yg[ng], sg[ng], ag, bg # Gaussian parameters + +int i, nx, xa, xb +real a, b, c, d, rad, wid, ya, yb, chisq +pointer xfit + +errchk id_mr_dofit + +begin + # Check starting values. + do i = 1, ng + if (IS_INDEF(x[i]) || (x[i] < 1) || (x[i] > npts)) + call error (1, "Invalid starting values") + + # Set minimum width and error radius. The minimum in the error radius + # is for defining the data window. The user error radius is used to + # check for an error in the derived center at the end of the centering. + + wid = max (width, MIN_WIDTH) + rad = max (2., radius) + + # Determine the pixel value range around the initial center, including + # the width and error radius buffer. Check for a minimum range. + + call alimr (x, ng, c, d) + x1 = max (1., c - wid / 2 - rad - wid) + x2 = min (real (npts), d + wid / 2 + rad + wid + 1) + nx = x2 - x1 + 1 + call alimr (data[x1], nx, a, b) + if (b - a < threshold) + call error (1, "Data range below threshold") + + # Allocate memory for the continuum subtracted data vector. The X + # range is just large enough to include the error radius and the + # half width. + + x1 = max (1., c - wid / 2 - rad) + x2 = min (real (npts), d + wid / 2 + rad + 1) + nx = x2 - x1 + 1 + + # Make the centering data positive, subtract the continuum, and + # apply a threshold to eliminate noise spikes. + + xa = nint(c) + ya = data[xa] + xb = nint(d) + yb = data[xb] + switch (type) { + case GEMISSION: + for (i = xa; i >= x1; i=i-1) + if (data[i] < ya) { + xa = i + ya = data[i] + } + for (i = xb; i <= x2; i=i+1) + if (data[i] < yb) { + xb = i + yb = data[i] + } + case GABSORPTION: + for (i = xa; i >= x1; i=i-1) + if (data[i] > ya) { + xa = i + ya = data[i] + } + for (i = xb; i <= x2; i=i+1) + if (data[i] > yb) { + xb = i + yb = data[i] + } + default: + call error (0, "Unknown feature type") + } + + # Set initial gaussian parameters. + if (IS_INDEF(ag) || IS_INDEF(bg)) { + if (xa == xb) + call error (1, "Can't determine background") + bg = (yb-ya) / (xb-xa) + ag = ya - bg * xa + } + do i = 1, ng { + xg[i] = x[i] + yg[i] = data[nint(x[i])] - ag - bg * x[i] + sg[i] = width / 6. + } + + # Determine the center. + call malloc (xfit, nx, TY_REAL) + do i = x1, x2 + Memr[xfit+i-x1] = i + + call id_mr_dofit (0, 3, 3, Memr[xfit], data[x1], nx, + ag, bg, xg, yg, sg, ng, chisq) + + # Check user centering error radius. + do i = 1, ng { + if (!IS_INDEF(xg[i])) { + if (abs (x[i] - xg[i]) > radius) + call error (2, "Error radius exceeded") + } + } + + # Free memory and return the center position. + call mfree (xfit, TY_REAL) +end diff --git a/noao/rv/rvidlines/idcolon.x b/noao/rv/rvidlines/idcolon.x new file mode 100644 index 00000000..c976c66b --- /dev/null +++ b/noao/rv/rvidlines/idcolon.x @@ -0,0 +1,288 @@ +include <gset.h> +include <error.h> +include <smw.h> +include "identify.h" + +# List of colon commands. +define CMDS "|show|features|image|nsum|database|read|write|add|coordlist|match\ + |maxfeatures|minsep|zwidth|labels|fwidth|ftype|cradius|threshold|" + +define SHOW 1 # Show parameters +define FEATURES 2 # Show list of features +define IMAGE 3 # Set new image +define NSUM 4 # Set the number of lines or columns to sum +define DATABASE 5 # Set new database +define READ 6 # Read database entry +define WRITE 7 # Write database entry +define ADD 8 # Add features from database +define COORDLIST 9 # Set new coordinate list +define MATCH 10 # Set coordinate list matching distance +define MAXFEATURES 11 # Set maximum number of features for auto find +define MINSEP 12 # Set minimum separation distance +define ZWIDTH 13 # Set zoom window width +define LABEL 14 # Set label type +define WIDTH 15 # Set centering width +define TYPE 16 # Set centering type +define RADIUS 17 # Set centering radius +define THRESHOLD 18 # Set the centering threshold + +# ID_COLON -- Respond to colon command. + +procedure id_colon (id, cmdstr, newimage, prfeature) + +pointer id # ID pointer +char cmdstr[ARB] # Colon command +char newimage[ARB] # New image name +int prfeature # Print current feature on status line + +char cmd[SZ_LINE] +int i, ncmd, ival[2] +real rval[2] +pointer im + +int nscan(), strdic() +pointer immap() +errchk immap, id_dbread, id_dbwrite, id_log + +begin + # Scan the command string and get the first word. + call sscan (cmdstr) + call gargwrd (cmd, SZ_LINE) + ncmd = strdic (cmd, cmd, SZ_LINE, CMDS) + + switch (ncmd) { + case SHOW: # :show - show values of parameters + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call gdeactivate (ID_GP(id), AW_CLEAR) + call id_show (id, "STDOUT") + call greactivate (ID_GP(id), AW_PAUSE) + } else { + iferr (call id_show (id, cmd)) { + call erract (EA_WARN) + prfeature = NO + } + } + case FEATURES: # :features - list features + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call gdeactivate (ID_GP(id), AW_CLEAR) + call id_log (id, "STDOUT", NULL) + call greactivate (ID_GP(id), AW_PAUSE) + } else { + iferr (call id_log (id, cmd, NULL)) { + call erract (EA_WARN) + prfeature = NO + } + } + case IMAGE: # :image - set image to identify + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call printf ("image %s\n") + call pargstr (Memc[ID_IMAGE(id)]) + prfeature = NO + } else { + call strcpy (cmd, newimage, SZ_FNAME) + iferr { + im = immap (newimage, READ_ONLY, 0) + call imunmap (im) + } then { + newimage[1] = EOS + call erract (EA_WARN) + prfeature = NO + } + } + case NSUM: # :nsum - set number of lines or columns to sum in image + call gargi (ival[1]) + if (nscan() == 1) { + call printf ("nsum %d %d\n") + call pargi (ID_NSUM(id,1)) + call pargi (ID_NSUM(id,2)) + prfeature = NO + } else { + ID_NSUM(id,1) = ival[1] + call gargi (ival[2]) + if (nscan() == 3) + ID_NSUM(id,2) = ival[2] + call smw_daxis (NULL, NULL, SMW_PAXIS(MW(ID_SH(id)),1), + ID_NSUM(id,1), ID_NSUM(id,2)) + } + case DATABASE: # :database - set database + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call printf ("database %s\n") + call pargstr (Memc[ID_DATABASE(id)]) + prfeature = NO + } else { + call strcpy (cmd, Memc[ID_DATABASE(id)], SZ_FNAME) + ID_NEWDBENTRY(id) = YES + } + case READ: # :read - read database entry + prfeature = NO + iferr { + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) + call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), + NO, YES) + else { + call gargi (ival[1]) + if (nscan() < 3) + ival[1] = ID_AP(id,1) + call gargi (ival[2]) + if (nscan() < 4) + ival[2] = ID_AP(id,2) + call id_dbread (id, cmd, ival, NO, YES) + } + } then + call erract (EA_WARN) + case WRITE: # :write - write database entry + prfeature = NO + iferr { + ival[1] = ID_AP(id,1) + ival[2] = ID_AP(id,2) + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) + call id_dbwrite (id, Memc[ID_IMAGE(id)], ival, YES) + else { + call gargi (ival[1]) + if (nscan() < 3) + ival[1] = ID_AP(id,1) + call gargi (ival[2]) + if (nscan() < 4) + ival[2] = ID_AP(id,2) + call id_dbwrite (id, cmd, ival, YES) + } + } then + call erract (EA_WARN) + case ADD: # :add - add features from database entry + prfeature = NO + iferr { + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) + call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), + YES, YES) + else { + call gargi (ival[1]) + if (nscan() < 3) + ival[1] = ID_AP(id,1) + call gargi (ival[2]) + if (nscan() < 4) + ival[2] = ID_AP(id,2) + call id_dbread (id, cmd, ival, YES, YES) + } + } then + call erract (EA_WARN) + case COORDLIST: # :coordlist - set coordinate list + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + call printf ("coordlist %s\n") + call pargstr (Memc[ID_COORDLIST(id)]) + prfeature = NO + } else { + call strcpy (cmd, Memc[ID_COORDLIST(id)], SZ_FNAME) + call id_unmapll (id) + call id_mapll (id) + } + case MATCH: # :match - set matching distance for coordinate list + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("match %g\n") + call pargr (ID_MATCH(id)) + prfeature = NO + } else + ID_MATCH(id) = rval[1] + case MAXFEATURES: # :maxfeatures - set max num features for auto find + call gargi (ival[1]) + if (nscan() == 1) { + call printf ("maxfeatures %d\n") + call pargi (ID_MAXFEATURES(id)) + prfeature = NO + } else + ID_MAXFEATURES(id) = ival[1] + case MINSEP: # :minsep - set minimum feature separation allowed + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("minsep %g\n") + call pargr (ID_MINSEP(id)) + prfeature = NO + } else + ID_MINSEP(id) = rval[1] + case ZWIDTH: # :zwidth - set zoom window width + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("zwidth %g\n") + call pargr (ID_ZWIDTH(id)) + prfeature = NO + } else { + ID_ZWIDTH(id) = rval[1] + if (ID_GTYPE(id) == 2) + ID_NEWGRAPH(id) = YES + } + case LABEL: # :labels - set label type + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + switch (ID_LABELS(id)) { + case 2: + call printf ("labels index\n") + case 3: + call printf ("labels pixel\n") + case 4: + call printf ("labels coords\n") + case 5: + call printf ("labels user\n") + case 6: + call printf ("labels both\n") + default: + call printf ("labels none\n") + } + prfeature = NO + } else { + ID_LABELS(id) = strdic (cmd, cmd, SZ_LINE, LABELS) + do i = 1, ID_NFEATURES(id) + call id_mark (id, i) + } + case WIDTH: # :fwidth - set centering width + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("fwidth %g\n") + call pargr (ID_FWIDTH(id)) + prfeature = NO + } else + ID_FWIDTH(id) = rval[1] + case TYPE: # :ftype - set centering type + call gargwrd (cmd, SZ_LINE) + if (nscan() == 1) { + switch (ID_FTYPE(id)) { + case EMISSION: + call printf ("ftype emission\n") + case ABSORPTION: + call printf ("ftype absorption\n") + case GEMISSION: + call printf ("ftype gemission\n") + case GABSORPTION: + call printf ("ftype gabsorption\n") + } + prfeature = NO + } else + ID_FTYPE(id) = strdic (cmd, cmd, SZ_LINE, FTYPES) + case RADIUS: # :cradius - set centering radius + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("cradius %g\n") + call pargr (ID_CRADIUS(id)) + prfeature = NO + } else + ID_CRADIUS(id) = rval[1] + case THRESHOLD: # :threshold - set centering threshold + call gargr (rval[1]) + if (nscan() == 1) { + call printf ("threshold %g\n") + call pargr (ID_THRESHOLD(id)) + prfeature = NO + } else + ID_THRESHOLD(id) = rval[1] + default: + call printf ("Unrecognized or ambiguous command\007") + prfeature = NO + } +end diff --git a/noao/rv/rvidlines/iddb.x b/noao/rv/rvidlines/iddb.x new file mode 100644 index 00000000..90bafea9 --- /dev/null +++ b/noao/rv/rvidlines/iddb.x @@ -0,0 +1,436 @@ +include <imset.h> +include <math/curfit.h> +include <smw.h> +include "identify.h" + +# ID_DBREAD -- Read features data from the database. + +procedure id_dbread (id, name, ap, add, verbose) + +pointer id # ID pointer +char name[SZ_LINE] # Image name +int ap[2] # Aperture number +int add # Add features? +int verbose # Verbose flag + +double pix +int i, j, k, ncoeffs, rec +pointer dt, sp, coeffs, line, str, sh + +int dtgeti(), dcvstati(), dtlocate(), dtscan(), nscan() +real dtgetr() +double dcvstatd() + +errchk dtremap(), dtlocate(), dtgeti(), dtgad() + +begin + call smark (sp) + call salloc (line, SZ_LINE, TY_CHAR) + + call strcpy ("id", Memc[line], SZ_LINE) + call imgcluster (name, Memc[line+2], SZ_LINE) + call dtremap (ID_DT(id), Memc[ID_DATABASE(id)], Memc[line], READ_ONLY) + + call id_dbsection (id, name, ap, Memc[ID_SECTION(id)], SZ_FNAME) + call sprintf (Memc[line], SZ_LINE, "identify %s%s") + call pargstr (name) + call pargstr (Memc[ID_SECTION(id)]) + + dt = ID_DT(id) + sh = ID_SH(id) + iferr (rec = dtlocate (dt, Memc[line])) { + call salloc (str, SZ_LINE, TY_CHAR) + call sprintf (Memc[str], SZ_LINE, "Entry not found: %s") + call pargstr (Memc[line]) + call error (0, Memc[str]) + } + + if (add == YES) { + j = dtgeti (dt, rec, "features") + k = j + ID_NFEATURES(id) + + call realloc (ID_PIX(id), k, TY_DOUBLE) + call realloc (ID_FIT(id), k, TY_DOUBLE) + call realloc (ID_USER(id), k, TY_DOUBLE) + call realloc (ID_WTS(id), k, TY_DOUBLE) + call realloc (ID_FWIDTHS(id), k, TY_REAL) + call realloc (ID_FTYPES(id), k, TY_INT) + call realloc (ID_LABEL(id), k, TY_POINTER) + + do i = 1, j { + k = dtscan (dt) + call gargd (pix) + + ID_NFEATURES(id) = ID_NFEATURES(id) + 1 + for (k=ID_NFEATURES(id); (k>1)&&(pix<PIX(id,k-1)); k=k-1) { + PIX(id,k) = PIX(id,k-1) + FIT(id,k) = FIT(id,k-1) + USER(id,k) = USER(id,k-1) + WTS(id,k) = WTS(id,k-1) + FWIDTH(id,k) = FWIDTH(id,k-1) + FTYPE(id,k) = FTYPE(id,k-1) + Memi[ID_LABEL(id)+k-1] = Memi[ID_LABEL(id)+k-2] + } + PIX(id,k) = pix + call gargd (FIT(id,k)) + call gargd (USER(id,k)) + call gargr (FWIDTH(id,k)) + call gargi (FTYPE(id,k)) + call gargd (WTS(id,k)) + call gargstr (Memc[line], SZ_LINE) + Memi[ID_LABEL(id)+k-1] = NULL + call id_label (Memc[line], Memi[ID_LABEL(id)+k-1]) + + # The following initialization is for backwards compatibility. + if (nscan() < 5) { + FWIDTH(id,k) = ID_FWIDTH(id) + FTYPE(id,k) = ID_FTYPE(id) + } else if (nscan() < 6) + WTS(id,k) = 1. + } + + } else { + if (SMW_FORMAT(MW(sh)) == SMW_ES || SMW_FORMAT(MW(sh)) == SMW_MS) { + iferr (APLOW(sh,1) = dtgetr (dt, rec, "aplow")) + APLOW(sh,1) = INDEF + iferr (APHIGH(sh,1) = dtgetr (dt, rec, "aphigh")) + APHIGH(sh,1) = INDEF + } + + do i = 1, ID_NFEATURES(id) + call mfree (Memi[ID_LABEL(id)+i-1], TY_CHAR) + + k = dtgeti (dt, rec, "features") + ID_NFEATURES(id) = k + ID_NALLOC(id) = k + call realloc (ID_PIX(id), k, TY_DOUBLE) + call realloc (ID_FIT(id), k, TY_DOUBLE) + call realloc (ID_USER(id), k, TY_DOUBLE) + call realloc (ID_WTS(id), k, TY_DOUBLE) + call realloc (ID_FWIDTHS(id), k, TY_REAL) + call realloc (ID_FTYPES(id), k, TY_INT) + call realloc (ID_LABEL(id), k, TY_POINTER) + + do i = 1, ID_NFEATURES(id) { + k = dtscan (dt) + call gargd (PIX(id,i)) + call gargd (FIT(id,i)) + call gargd (USER(id,i)) + call gargr (FWIDTH(id,i)) + call gargi (FTYPE(id,i)) + call gargd (WTS(id,i)) + call gargstr (Memc[line], SZ_LINE) + Memi[ID_LABEL(id)+i-1] = NULL + call id_label (Memc[line], Memi[ID_LABEL(id)+i-1]) + + # The following initialization is for backwards compatibility. + if (nscan() < 5) { + FWIDTH(id,i) = ID_FWIDTH(id) + FTYPE(id,i) = ID_FTYPE(id) + } else if (nscan() < 6) + WTS(id,i) = 1. + } + + iferr (ID_SHIFT(id) = dtgetr (dt, rec, "shift")) + ID_SHIFT(id) = 0. + iferr (ID_REDSHIFT(id) = dtgetr (dt, rec, "redshift")) + ID_REDSHIFT(id) = 0. + iferr (ID_RMSRED(id) = dtgetr (dt, rec, "redshift_rms")) + ID_RMSRED(id) = 0. + + iferr { + ncoeffs = dtgeti (dt, rec, "coefficients") + call salloc (coeffs, ncoeffs, TY_DOUBLE) + call dtgad (dt, rec, "coefficients", Memd[coeffs], ncoeffs, + ncoeffs) + + if (ID_CV(id) != NULL) + call dcvfree (ID_CV(id)) + call dcvrestore (ID_CV(id), Memd[coeffs]) + + call ic_putr (ID_IC(id), "xmin", real (dcvstatd(ID_CV(id), + CVXMIN))) + call ic_putr (ID_IC(id), "xmax", real (dcvstatd(ID_CV(id), + CVXMAX))) + ifnoerr (call dtgstr (dt,rec,"function",Memc[line],SZ_LINE)) { + call ic_pstr (ID_IC(id), "function", Memc[line]) + call ic_puti (ID_IC(id), "order", dtgeti (dt, rec, "order")) + call dtgstr (dt, rec, "sample", Memc[line], SZ_LINE) + call ic_pstr (ID_IC(id), "sample", Memc[line]) + call ic_puti (ID_IC(id), "naverage", + dtgeti (dt, rec, "naverage")) + call ic_puti (ID_IC(id), "niterate", + dtgeti (dt, rec, "niterate")) + call ic_putr (ID_IC(id), "low", + dtgetr (dt, rec, "low_reject")) + call ic_putr (ID_IC(id), "high", + dtgetr (dt, rec, "high_reject")) + call ic_putr (ID_IC(id), "grow", dtgetr (dt, rec, "grow")) + } else { + call ic_puti (ID_IC(id), "order", dcvstati (ID_CV(id), + CVORDER)) + switch (dcvstati (ID_CV(id), CVTYPE)) { + case LEGENDRE: + call ic_pstr (ID_IC(id), "function", "legendre") + case CHEBYSHEV: + call ic_pstr (ID_IC(id), "function", "chebyshev") + case SPLINE1: + call ic_pstr (ID_IC(id), "function", "spline1") + case SPLINE3: + call ic_pstr (ID_IC(id), "function", "spline3") + } + } + + ID_NEWCV(id) = YES + ID_CURRENT(id) = min (1, ID_NFEATURES(id)) + } then + ; + } + + call sfree (sp) + + if (ID_NFEATURES(id) > 0) { + ID_NEWGRAPH(id) = YES + ID_NEWFEATURES(id) = YES + ID_CURRENT(id) = 1 + } else + ID_CURRENT(id) = 0 + + if (verbose == YES) { + call printf ("identify %s%s\n") + call pargstr (name) + call pargstr (Memc[ID_SECTION(id)]) + } +end + + +# ID_DBWRITE -- Write features data to the database. + +procedure id_dbwrite (id, name, ap, verbose) + +pointer id # ID pointer +char name[ARB] # Image name +int ap[2] # Aperture number +int verbose # Verbose flag + +int i, ncoeffs +pointer dt, sp, coeffs, root, sh, im + +int dcvstati(), ic_geti() +real ic_getr() + +errchk dtremap + +begin + call smark (sp) + call salloc (root, SZ_FNAME, TY_CHAR) + + call strcpy ("id", Memc[root], SZ_FNAME) + call imgcluster (name, Memc[root+2], SZ_FNAME) + call dtremap (ID_DT(id), Memc[ID_DATABASE(id)], Memc[root], APPEND) + + call id_dbsection (id, name, ap, Memc[ID_SECTION(id)], SZ_FNAME) + + sh = ID_SH(id) + dt = ID_DT(id) + call dtptime (dt) + call dtput (dt, "begin\tidentify %s%s\n") + call pargstr (name) + call pargstr (Memc[ID_SECTION(id)]) + call dtput (dt, "\tid\t%s\n") + call pargstr (name) + call dtput (dt, "\ttask\tidentify\n") + call dtput (dt, "\timage\t%s%s\n") + call pargstr (Memc[ID_IMAGE(id)]) + call pargstr (Memc[ID_SECTION(id)]) + if (SMW_FORMAT(MW(sh)) == SMW_ES || SMW_FORMAT(MW(sh)) == SMW_MS) { + call dtput (dt, "\taperture\t%d\n") + call pargi (ID_AP(id,1)) + call dtput (dt, "\taplow\t%g\n") + call pargr (APLOW(sh,1)) + call dtput (dt, "\taphigh\t%g\n") + call pargr (APHIGH(sh,1)) + } + + call dtput (dt, "\tfeatures\t%d\n") + call pargi (ID_NFEATURES(id)) + do i = 1, ID_NFEATURES(id) { + call dtput (dt, "\t %10.2f %10.8g %10.8g %5.1f %d %d %s\n") + call pargd (PIX(id,i)) + call pargd (FIT(id,i)) + call pargd (USER(id,i)) + call pargr (FWIDTH(id,i)) + call pargi (FTYPE(id,i)) + call pargd (WTS(id,i)) + if (Memi[ID_LABEL(id)+i-1] != NULL) + call pargstr (Memc[Memi[ID_LABEL(id)+i-1]]) + else + call pargstr ("") + } + + if (ID_SHIFT(id) != 0.) { + call dtput (dt, "\tshift\t%g\n") + call pargd (ID_SHIFT(id)) + } + if (ID_REDSHIFT(id) != 0.) { + call dtput (dt, "\tredshift\t%g\n") + call pargd (ID_REDSHIFT(id)) + call dtput (dt, "\tredshift_rms\t%g\n") + call pargd (ID_RMSRED(id)) + } + + if (ID_CV(id) != NULL) { + call dtput (dt, "\tfunction %s\n") + call ic_gstr (ID_IC(id), "function", Memc[root], SZ_FNAME) + call pargstr (Memc[root]) + call dtput (dt, "\torder %d\n") + call pargi (ic_geti (ID_IC(id), "order")) + call dtput (dt, "\tsample %s\n") + call ic_gstr (ID_IC(id), "sample", Memc[root], SZ_FNAME) + call pargstr (Memc[root]) + call dtput (dt, "\tnaverage %d\n") + call pargi (ic_geti (ID_IC(id), "naverage")) + call dtput (dt, "\tniterate %d\n") + call pargi (ic_geti (ID_IC(id), "niterate")) + call dtput (dt, "\tlow_reject %g\n") + call pargr (ic_getr (ID_IC(id), "low")) + call dtput (dt, "\thigh_reject %g\n") + call pargr (ic_getr (ID_IC(id), "high")) + call dtput (dt, "\tgrow %g\n") + call pargr (ic_getr (ID_IC(id), "grow")) + + ncoeffs = dcvstati (ID_CV(id), CVNSAVE) + call salloc (coeffs, ncoeffs, TY_DOUBLE) + call dcvsave (ID_CV(id), Memd[coeffs]) + call dtput (dt, "\tcoefficients\t%d\n") + call pargi (ncoeffs) + do i = 1, ncoeffs { + call dtput (dt, "\t\t%g\n") + call pargd (Memd[coeffs+i-1]) + } + } + + call dtput (dt, "\n") + + ID_NEWFEATURES(id) = NO + ID_NEWCV(id) = NO + ID_NEWDBENTRY(id) = NO + + if (ID_TASK(id) == IDENTIFY) { + if (verbose == YES) { + call printf ("identify %s%s\n") + call pargstr (name) + call pargstr (Memc[ID_SECTION(id)]) + } + + # Enter reference spectrum name in image header. + im = IM(sh) + call imseti (im, IM_WHEADER, YES) + call imastr (im, "REFSPEC1", Memc[ID_IMAGE(id)]) + iferr (call imdelf (im, "REFSPEC2")) + ; + } + + call sfree (sp) +end + + +# ID_DBCHECK -- Check if there is an entry in the database. +# This does not actually read the database entry. It also assumes that +# if a database is already open it is for the same image (the image +# names are not checked) and the database has been scanned. + +int procedure id_dbcheck (id, name, ap) + +pointer id # ID pointer +char name[SZ_LINE] # Image name +int ap[2] # Aperture number + +int rec, stat +pointer sp, line, sec + +int dtlocate() + +errchk dtremap(), dtlocate() + +begin + call smark (sp) + call salloc (line, SZ_LINE, TY_CHAR) + call salloc (sec, SZ_LINE, TY_CHAR) + + if (ID_DT(id) == NULL) { + call strcpy ("id", Memc[line], SZ_LINE) + call imgcluster (name, Memc[line+2], SZ_LINE) + iferr (call dtremap (ID_DT(id), Memc[ID_DATABASE(id)], Memc[line], + READ_ONLY)) { + call sfree (sp) + return (NO) + } + } + + call id_dbsection (id, name, ap, Memc[sec], SZ_LINE) + call sprintf (Memc[line], SZ_LINE, "identify %s%s") + call pargstr (name) + call pargstr (Memc[sec]) + + iferr (rec = dtlocate (ID_DT(id), Memc[line])) + stat = NO + else + stat = YES + + call sfree (sp) + return (stat) +end + + +# ID_DBSECTION -- Make the IDENTIFY section. + +procedure id_dbsection (id, name, ap, section, sz_section) + +pointer id #I ID pointer +char name[SZ_LINE] #I Image name +int ap[2] #I Aperture number +char section[sz_section] #O IDENTIFY section +int sz_section #I Size of section string + +pointer sh, smw +bool streq() + +begin + sh = ID_SH(id) + smw = MW(sh) + + switch (SMW_FORMAT(smw)) { + case SMW_ND: + section[1] = EOS + if (streq (name, Memc[ID_IMAGE(id)])) { + switch (SMW_LDIM(smw)) { + case 2: + switch (SMW_LAXIS(smw,1)) { + case 1: + call sprintf (section, sz_section, "[*,%d]") + case 2: + call sprintf (section, sz_section, "[%d,*]") + } + #call pargi (LINDEX(sh,1)) + call pargi (ap[1]) + case 3: + switch (SMW_LAXIS(smw,1)) { + case 1: + call sprintf (section, sz_section, "[*,%d,%d]") + case 2: + call sprintf (section, sz_section, "[%d,*,%d]") + case 3: + call sprintf (section, sz_section, "[%d,%d,*]") + } + #call pargi (LINDEX(sh,1)) + #call pargi (LINDEX(sh,2)) + call pargi (ap[1]) + call pargi (ap[2]) + } + } + case SMW_ES, SMW_MS: + call sprintf (section, sz_section, " - Ap %d") + call pargi (ap[1]) + } +end diff --git a/noao/rv/rvidlines/iddeblend.x b/noao/rv/rvidlines/iddeblend.x new file mode 100644 index 00000000..75e32a78 --- /dev/null +++ b/noao/rv/rvidlines/iddeblend.x @@ -0,0 +1,413 @@ +include <mach.h> + + +# ID_MR_DOFIT -- Fit gaussian components. This is an interface to ID_DOFIT1 +# which puts parameters into the form required by ID_DOFIT1 and vice-versa. +# It also implements a constrained approach to the solution. + +procedure id_mr_dofit (bkgfit, posfit, sigfit, x, y, npts, y1, dy, xg, yg, sg, + ng, chisq) + +int bkgfit # Fit background (0=no, 1=yes) +int posfit # Position fitting flag (1=fixed, 2=single, 3=all) +int sigfit # Sigma fitting flag (1=fixed, 2=single, 3=all) +real x[npts] # X data +real y[npts] # Y data +int npts # Number of points +real y1 # Continuum offset +real dy # Continuum slope +real xg[ng] # Initial and final x coordinates of gaussians +real yg[ng] # Initial and final y coordinates of gaussians +real sg[ng] # Initial and final sigmas of gaussians +int ng # Number of gaussians +real chisq # Chi squared + +int i +pointer sp, a, j +errchk id_dofit1 + +begin + call smark (sp) + call salloc (a, 4 + 3 * ng, TY_REAL) + + # Convert positions and widths relative to first component. + Memr[a] = y1 + Memr[a+1] = dy + Memr[a+2] = xg[1] + Memr[a+3] = sg[1] + do i = 1, ng { + j = a + 3 * i + 1 + Memr[j] = yg[i] + Memr[j+1] = xg[i] - Memr[a+2] + Memr[j+2] = sg[i] / Memr[a+3] + } + + # Do fit. + do i = 0, bkgfit { + switch (10*posfit+sigfit) { + case 11: + call id_dofit1 (i, 1, 1, x, y, npts, Memr[a], ng, chisq) + case 12: + call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq) + case 13: + call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 1, 3, x, y, npts, Memr[a], ng, chisq) + case 21: + call id_dofit1 (i, 2, 1, x, y, npts, Memr[a], ng, chisq) + case 22: + call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 2, 2, x, y, npts, Memr[a], ng, chisq) + case 23: + call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 2, 2, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 2, 3, x, y, npts, Memr[a], ng, chisq) + case 31: + call id_dofit1 (i, 2, 1, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 3, 1, x, y, npts, Memr[a], ng, chisq) + case 32: + call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 2, 2, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 3, 2, x, y, npts, Memr[a], ng, chisq) + case 33: + call id_dofit1 (i, 1, 2, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 2, 2, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 3, 2, x, y, npts, Memr[a], ng, chisq) + call id_dofit1 (i, 3, 3, x, y, npts, Memr[a], ng, chisq) + } + } + + y1 = Memr[a] + dy = Memr[a+1] + do i = 1, ng { + j = a + 3 * i + 1 + yg[i] = Memr[j] + xg[i] = Memr[j+1] + Memr[a+2] + sg[i] = abs (Memr[j+2] * Memr[a+3]) + } + + call sfree (sp) +end + + +# ID_MODEL -- Compute model. +# +# I(x) = I(i) exp {-[(x-xg(i)) / sg(i)]**2 / 2.} +# +# where the params are I1, I2, xg, yg, and sg. + +real procedure id_model (x, xg, yg, sg, ng) + +real x # X value to be evaluated +real xg[ng] # X coordinates of gaussians +real yg[ng] # Y coordinates of gaussians +real sg[ng] # Sigmas of gaussians +int ng # Number of gaussians + +int i +real y, arg + +begin + y = 0. + do i = 1, ng { + arg = (x - xg[i]) / sg[i] + if (abs (arg) < 7.) + y = y + yg[i] * exp (-arg**2 / 2.) + } + return (y) +end + + +# ID_DOFIT1 -- Perform nonlinear iterative fit for the specified parameters. +# This uses the Levenberg-Marquardt method from NUMERICAL RECIPES. + +procedure id_dofit1 (bkgfit, posfit, sigfit, x, y, npts, a, nlines, chisq) + +int bkgfit # Background fit (0=no, 1=yes) +int posfit # Position fitting flag (1=fixed, 2=one, 3=all) +int sigfit # Sigma fitting flag (1=fixed, 2=one, 3=all) +real x[npts] # X data +real y[npts] # Y data +int npts # Number of points +real a[ARB] # Fitting parameters +int nlines # Number of lines +real chisq # Chi squared + +int i, np, nfit +real mr, chi2 +pointer sp, flags, ptr +errchk id_mr_solve + +begin + # Number of terms is 3 for each line plus common background, center + # and sigma. + + np = 3 * nlines + 4 + + call smark (sp) + call salloc (flags, np, TY_INT) + ptr = flags + + # Background. + if (bkgfit == 1) { + Memi[ptr] = 1 + Memi[ptr+1] = 2 + ptr = ptr + 2 + } + + # Peaks are always fit. + do i = 1, nlines { + Memi[ptr] = 3 * i + 2 + ptr = ptr + 1 + } + + # Positions. + switch (posfit) { + case 2: + Memi[ptr] = 3 + ptr = ptr + 1 + case 3: + Memi[ptr] = 3 + ptr = ptr + 1 + do i = 1, nlines { + Memi[ptr] = 3 * i + 3 + ptr = ptr + 1 + } + } + + # Sigmas. + switch (sigfit) { + case 2: + Memi[ptr] = 4 + ptr = ptr + 1 + case 3: + Memi[ptr] = 4 + ptr = ptr + 1 + do i = 1, nlines { + Memi[ptr] = 3 * i + 4 + ptr = ptr + 1 + } + } + + nfit = ptr - flags + mr = -1. + i = 0 + chi2 = MAX_REAL + repeat { + call id_mr_solve (x, y, npts, a, Memi[flags], np, nfit, mr, chisq) + if (chi2 - chisq > 1.) + i = 0 + else + i = i + 1 + chi2 = chisq + } until (i == 3) + + mr = 0. + call id_mr_solve (x, y, npts, a, Memi[flags], np, nfit, mr, chisq) + + call sfree (sp) +end + + +# ID_DERIVS -- Compute model and derivatives for MR_SOLVE procedure. +# +# I(x) = I1 + I2 * x + I(i) exp {-[(x-xc-dx(i)) / (sig * sig(i))]**2 / 2.} +# +# where the params are I1, I2, xc, sig, I(i), dx(i), and sig(i) (i=1,nlines). + +procedure id_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 + +int i +real sig, arg, ex, fac + +begin + y = a[1] + a[2] * x + dyda[1] = 1. + dyda[2] = x + dyda[3] = 0. + dyda[4] = 0. + do i = 5, na, 3 { + sig = a[4] * a[i+2] + arg = (x - a[3] - a[i+1]) / sig + if (abs (arg) < 7.) + ex = exp (-arg**2 / 2.) + else + ex = 0. + fac = a[i] * ex * arg + + y = y + a[i] * ex + dyda[3] = dyda[3] + fac / sig + dyda[4] = dyda[4] + fac * arg / a[4] + dyda[i] = ex + dyda[i+1] = fac / sig + dyda[i+2] = fac * arg / a[i+2] + } +end + + +# ID_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 id_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 id_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 id_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 id_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 id_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 + + +# ID_MR_EVAL -- Evaluate curvature matrix. This calls procedure DERIVS. + +procedure id_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 id_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 id_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/noao/rv/rvidlines/iddelete.x b/noao/rv/rvidlines/iddelete.x new file mode 100644 index 00000000..cd96abb1 --- /dev/null +++ b/noao/rv/rvidlines/iddelete.x @@ -0,0 +1,26 @@ +include "identify.h" + +# ID_DELETE -- Delete a feature. + +procedure id_delete (id, feature) + +pointer id # ID pointer +int feature # Feature to be deleted + +int i + +begin + call mfree (Memi[ID_LABEL(id)+feature-1], TY_CHAR) + do i = feature + 1, ID_NFEATURES(id) { + PIX(id,i-1) = PIX(id,i) + FIT(id,i-1) = FIT(id,i) + USER(id,i-1) = USER(id,i) + WTS(id,i-1) = WTS(id,i) + FWIDTH(id,i-1) = FWIDTH(id,i) + FTYPE(id,i-1) = FTYPE(id,i) + Memi[ID_LABEL(id)+i-2] = Memi[ID_LABEL(id)+i-1] + } + Memi[ID_LABEL(id)+ID_NFEATURES(id)-1] = NULL + ID_NFEATURES(id) = ID_NFEATURES(id) - 1 + ID_NEWFEATURES(id) = YES +end diff --git a/noao/rv/rvidlines/iddofit.x b/noao/rv/rvidlines/iddofit.x new file mode 100644 index 00000000..691e7350 --- /dev/null +++ b/noao/rv/rvidlines/iddofit.x @@ -0,0 +1,101 @@ +include "identify.h" + +# ID_DOFIT -- Fit a function to the features. Eliminate INDEF points. + +procedure id_dofit (id, interactive) + +pointer id # ID pointer +int interactive # Interactive fit? + +int i, j, k, nfit, ic_geti() +pointer gt1, sp, x, y, wts, rejpts, str, gt_init() + +begin + if (ID_NFEATURES(id) == 0) { + if (ID_CV(id) != NULL) { + call dcvfree (ID_CV(id)) + ID_SHIFT(id) = 0. + ID_REDSHIFT(id) = 0. + ID_NEWGRAPH(id) = YES + ID_NEWCV(id) = YES + } + return + } + + call smark (sp) + call salloc (x, ID_NFEATURES(id), TY_DOUBLE) + call salloc (y, ID_NFEATURES(id), TY_DOUBLE) + call salloc (wts, ID_NFEATURES(id), TY_DOUBLE) + + nfit = 0 + do i = 1, ID_NFEATURES(id) { + if (IS_INDEFD (PIX(id,i)) || IS_INDEFD (USER(id,i))) + next + Memd[x+nfit] = PIX(id,i) + Memd[y+nfit] = USER(id,i) + Memd[wts+nfit] = 1. + nfit = nfit + 1 + } + + if (nfit > 1) { + if (interactive == YES) { + call salloc (str, SZ_LINE, TY_CHAR) + gt1 = gt_init() + call icg_fitd (ID_IC(id), ID_GP(id), "cursor", gt1, ID_CV(id), + Memd[x], Memd[y], Memd[wts], nfit) + call gt_free (gt1) + } else + call ic_fitd (ID_IC(id), ID_CV(id), Memd[x], Memd[y], Memd[wts], + nfit, YES, YES, YES, YES) + + if (ic_geti (ID_IC(id), "nreject") > 0 && + ic_geti (ID_IC(id), "nfit") == nfit) + rejpts = ic_geti (ID_IC(id), "rejpts") + else + rejpts = NULL + + j = 0 + k = 0 + do i = 1, ID_NFEATURES(id) { + if (IS_INDEFD (PIX(id,i)) || IS_INDEFD (USER(id,i))) { + j = j + 1 + PIX(id,j) = PIX(id,i) + FIT(id,j) = FIT(id,i) + USER(id,j) = USER(id,i) + WTS(id,j) = WTS(id,i) + FWIDTH(id,j) = FWIDTH(id,i) + FTYPE(id,j) = FTYPE(id,i) + } else { + if (Memd[wts+k] != 0.) { + j = j + 1 + PIX(id,j) = Memd[x+k] + FIT(id,j) = FIT(id,i) + USER(id,j) = Memd[y+k] + WTS(id,j) = Memd[wts+k] + if (rejpts != NULL) + if (Memi[rejpts+k] == YES) + WTS(id,j) = 0. + FWIDTH(id,j) = FWIDTH(id,i) + FTYPE(id,j) = FTYPE(id,i) + } + k = k + 1 + } + } + ID_NFEATURES(id) = j + + ID_SHIFT(id) = 0. + ID_REDSHIFT(id) = 0. + ID_NEWCV(id) = YES + ID_NEWGRAPH(id) = YES + } else { + if (ID_CV(id) != NULL) { + call dcvfree (ID_CV(id)) + ID_SHIFT(id) = 0. + ID_REDSHIFT(id) = 0. + ID_NEWCV(id) = YES + ID_NEWGRAPH(id) = YES + } + } + + call sfree (sp) +end diff --git a/noao/rv/rvidlines/iddoshift.x b/noao/rv/rvidlines/iddoshift.x new file mode 100644 index 00000000..6e397455 --- /dev/null +++ b/noao/rv/rvidlines/iddoshift.x @@ -0,0 +1,42 @@ +include "identify.h" + +# ID_DOSHIFT -- Minimize residuals by constant shift. + +procedure id_doshift (id, interactive) + +pointer id # ID pointer +int interactive # Called interactively? + +int i, j +double shft, delta, rms, id_fitpt() + +begin + shft = 0. + rms = 0. + j = 0 + for (i=1; i <= ID_NFEATURES(id); i = i + 1) { + if (IS_INDEFD (USER(id,i)) || WTS(id,i) == 0.) + next + delta = USER(id,i) - id_fitpt (id, PIX(id,i)) + shft = shft + delta + rms = rms + delta * delta + j = j + 1 + } + + if (j > 0) { + shft = shft / j + rms = rms / j + if (interactive == YES) { + call printf ("%s%s: Coordinate shift=%5f, rms=%5f, npts=%3d\n") + call pargstr (Memc[ID_IMAGE(id)]) + call pargstr (Memc[ID_SECTION(id)]) + call pargd (shft) + call pargd (sqrt (rms - shft ** 2)) + call pargi (j) + } + ID_SHIFT(id) = ID_SHIFT(id) + shft + ID_REDSHIFT(id) = 0. + ID_NEWCV(id) = YES + ID_NEWGRAPH(id) = YES + } +end diff --git a/noao/rv/rvidlines/identify.h b/noao/rv/rvidlines/identify.h new file mode 100644 index 00000000..ad34b767 --- /dev/null +++ b/noao/rv/rvidlines/identify.h @@ -0,0 +1,97 @@ +# Task parameters + +define LEN_IDSTRUCT 64 # Length ID structure + +define ID_TASK Memi[$1] # Task ID +define ID_IMAGE Memi[$1+1] # Image (pointer) +define ID_SECTION Memi[$1+2] # Section for 2D and 3D images (pointer) +define ID_LINE Memi[$1+$2+2] # Image line or column [2] +define ID_MAXLINE Memi[$1+$2+4] # Maximum line or column [2] +define ID_AP Memi[$1+$2+6] # Aperture if appropriate [2] +define ID_APS Memi[$1+9] # Array of apertures (pointer) +define ID_NSUM Memi[$1+$2+10] # Number of lines to sum [2] +define ID_MAXFEATURES Memi[$1+13] # Maximum number of features +define ID_FTYPE Memi[$1+14] # Feature type +define ID_MINSEP Memr[P2R($1+15)] # Minimum pixel separation +define ID_MATCH Memr[P2R($1+16)] # Maximum matching separation +define ID_FWIDTH Memr[P2R($1+17)] # Feature width in pixels +define ID_CRADIUS Memr[P2R($1+18)] # Centering radius in pixels +define ID_THRESHOLD Memr[P2R($1+19)] # Centering threshold +define ID_ZWIDTH Memr[P2R($1+20)] # Zoom window width in fit units +define ID_DATABASE Memi[$1+21] # Name of database (pointer) +define ID_COORDLIST Memi[$1+22] # Name of coordinate list (pointer) +define ID_LL Memi[$1+23] # Pointer to lines in coordinate list +define ID_LABELS Memi[$1+24] # Type of feature labels +define ID_LOGFILES Memi[$1+25] # List of logfiles + +# Common image data + +define ID_SHIFT Memd[P2D($1+26)]# Wavelength shift +define ID_REDSHIFT Memd[P2D($1+28)]# Redshift of spectrum +define ID_RMSRED Memd[P2D($1+30)]# Redshift of spectrum +define ID_ZHELIO Memd[P2D($1+32)]# Heliocentric correction in redshift +define ID_IMDATA Memi[$1+34] # Image data (pointer) +define ID_PIXDATA Memi[$1+35] # Pixel coordinates (pointer) +define ID_FITDATA Memi[$1+36] # Fit coordinates (pointer) +define ID_NPTS Memi[$1+37] # Number of points + +# Features + +define ID_NFEATURES Memi[$1+38] # Number of features +define ID_NALLOC Memi[$1+39] # Length of allocated feature arrays +define ID_PIX Memi[$1+40] # Feature pixel coordinates (pointer) +define ID_FIT Memi[$1+41] # Feature fit coordinates (pointer) +define ID_USER Memi[$1+42] # Feature user coordinates (pointer) +define ID_WTS Memi[$1+43] # Feature weights (pointer) +define ID_FWIDTHS Memi[$1+44] # Feature width (pointer) +define ID_FTYPES Memi[$1+45] # Feature type (pointer) +define ID_LABEL Memi[$1+46] # Feature label (pointer) +define ID_CURRENT Memi[$1+47] # Current feature + +# Pointers for other packages and to save data + +define ID_SH Memi[$1+48] # SHDR pointer +define ID_LP Memi[$1+49] # Logical to physical transformation +define ID_PL Memi[$1+50] # Physical to logical transformation +define ID_IC Memi[$1+51] # ICFIT pointer +define ID_CV Memi[$1+52] # Curfit pointer +define ID_GP Memi[$1+53] # GIO pointer +define ID_GT Memi[$1+54] # Gtools pointer +define ID_ID Memi[$1+55] # Array of structure pointers (pointer) +define ID_NID Memi[$1+56] # Number of saved structure +define ID_DT Memi[$1+57] # Database pointer + +# Flags + +define ID_NEWFEATURES Memi[$1+58] # Has feature list changed? +define ID_NEWCV Memi[$1+59] # Has fitting function changed? +define ID_NEWGRAPH Memi[$1+60] # Has graph changed? +define ID_NEWDBENTRY Memi[$1+61] # Has database entry changed? +define ID_REFIT Memi[$1+62] # Refit feature data? +define ID_GTYPE Memi[$1+63] # Graph type + +# End of structure ---------------------------------------------------------- + +# Task ID +define IDENTIFY 1 # Standard identify +define RVIDLINES 2 # Line radial velocities + +define LABELS "|none|index|pixel|coord|user|both|" +define FTYPES "|emission|absorption|gemission|gabsorption|" +define EMISSION 1 # Emission feature (center1d) +define ABSORPTION 2 # Absorption feature (center1d) +define GEMISSION 3 # Emission feature (center1d) +define GABSORPTION 4 # Absorption feature (center1d) + +define IMDATA Memr[ID_IMDATA($1)+$2-1] +define PIXDATA Memd[ID_PIXDATA($1)+$2-1] +define FITDATA Memd[ID_FITDATA($1)+$2-1] + +define PIX Memd[ID_PIX($1)+$2-1] +define FIT Memd[ID_FIT($1)+$2-1] +define USER Memd[ID_USER($1)+$2-1] +define WTS Memd[ID_WTS($1)+$2-1] +define FWIDTH Memr[ID_FWIDTHS($1)+$2-1] +define FTYPE Memi[ID_FTYPES($1)+$2-1] + +define VLIGHT 2.997925e5 # Speed of light, Km/sec diff --git a/noao/rv/rvidlines/identify.key b/noao/rv/rvidlines/identify.key new file mode 100644 index 00000000..ae2f9339 --- /dev/null +++ b/noao/rv/rvidlines/identify.key @@ -0,0 +1,104 @@ + IDENTIFY HELP + + +STATUS LINE + +The status line gives the pixel position, fitted wavelength, user wavelength, +wavelength residual, and optional line identification: + + pixel fitted user residual [identification] + + +CURSOR KEY SUMMARY + +? Help l Match list (fit) w Window graph +a Affect all features m Mark feature x Crosscorrelate peaks +b Deblend n Next feature y Find peaks +c Center feature(s) o Go to line z Zoom graph +d Delete feature(s) p Pan graph + Next feature +f Fit positions q Quit - Previous feature +g Fit zero point shift r Redraw graph . Nearest feature +i Initialize s Shift feature I Interrupt +j Preceding line t Reset position +k Next line u Enter coordinate + + +COLON COMMAND SUMMARY + +:show [file] :features [file] :coordlist [file] +:cradius [value] :threshold [value] :database [file] +:ftype [type] :fwidth [value] :image [image] +:labels [type] :match [value] :maxfeatures [value] +:minsep [value] :read [image [ap]] :write [image [ap]] +:zwidth [value] :add [image [ap]] :vlog [file] + + +CURSOR KEYS + +? Clear the screen and print menu of options +a Apply next (c)enter or (d)elete operation to (a)ll features +b Mark and de(b)lend features by Gaussian fitting +c (C)enter the feature nearest the cursor +d (D)elete the feature nearest the cursor +f (F)it a function of pixel coordinate to the user coordinates +g Fit a zero point shift to the user coordinates +h Match coordinates in the coordinate list without modifying the fit +i (I)nitialize (delete features and coordinate fit) +j Go to the preceding image line or column in a 2D or multispec image +k Go to the next image line or column in a 2D or multispec image +l Match coordinates in the coordinate (l)ist and fit/refit the dispersion +m (M)ark a new feature near the cursor and enter coordinate and label +n Move the cursor or zoom to the (n)ext feature (same as +) +o Go to the specified image line or column in a 2D or multispec image +p (P)an to user defined window after (z)ooming on a feature +q (Q)uit and continue with next image (also carriage return) +r (R)edraw the graph +s (S)hift the current feature to the position of the cursor +t Reset the position of a feature without centering +u Enter a new (u)ser coordinate and label for the current feature +v Compute a redshift and velocity from the fitted and user coordinates +w (W)indow the graph. Use '?' to window prompt for more help. +x Crosscorrelate features with the data peaks and reregister +y Automatically find "maxfeatures" strongest peaks and identify them +z (Z)oom on the feature nearest the cursor +. Move the cursor or zoom to the feature nearest the cursor ++ Move the cursor or zoom to the next feature +- Move the cursor or zoom to the previous feature +I Interrupt task and exit immediately. Database information is not saved. + + +COLON COMMANDS + +The parameters are listed or set with the following commands which may be +abbreviated. To list the value of a parameter type the command alone. + +:show file Show the values of all the parameters +:features file Write feature list to file (default is STDOUT) +:vlog file Write velocity information to file (default is STDOUT) + +:coordlist file Coordinate list file +:cradius value Centering radius in pixels +:threshold value Detection threshold for feature centering +:database name Database for recording feature records +:ftype value Feature type (emission or absorption) +:fwidth value Feature width in pixels +:image imagename Set a new image or show the current image +:labels value Feature label type (none|index|pixel|coords|user|both) +:match value Coordinate list matching distance +:maxfeatures value Maximum number of features automatically found +:minsep value Minimum separation allowed between features +:read name ap Read a record from the database + (name and ap default to the current spectrum) +:write name ap Write a record to the database + (name and ap default to the current spectrum) +:add name ap Add features from the database + (name and ap default to the current spectrum) +:zwidth value Zoom width in user units + +Labels: + none - No labels + index - Sequential numbers in order of increasing pixel position + pixel - Pixel coordinates + coords - User coordinates such as wavelength + user - User labels + both - Combination of coords and user diff --git a/noao/rv/rvidlines/idfitdata.x b/noao/rv/rvidlines/idfitdata.x new file mode 100644 index 00000000..48c36984 --- /dev/null +++ b/noao/rv/rvidlines/idfitdata.x @@ -0,0 +1,140 @@ +include <smw.h> +include "identify.h" + +# ID_FITDATA -- Compute fit coordinates from pixel coordinates. + +procedure id_fitdata (id) + +pointer id # ID pointer +int i + +begin + call mfree (ID_FITDATA(id), TY_DOUBLE) + call malloc (ID_FITDATA(id), ID_NPTS(id), TY_DOUBLE) + + if (ID_CV(id) == NULL) + call achtrd (Memr[SX(ID_SH(id))], FITDATA(id,1), ID_NPTS(id)) + else { + call dcvvector (ID_CV(id), PIXDATA(id,1), FITDATA(id,1), + ID_NPTS(id)) + if (FITDATA(id,2) > FITDATA(id,1)) { + do i = 3, ID_NPTS(id) + if (FITDATA(id,i) < FITDATA(id,i-1)) + call error (1, "Coordinate solution is not monotonic") + } else { + do i = 3, ID_NPTS(id) + if (FITDATA(id,i) > FITDATA(id,i-1)) + call error (1, "Coordinate solution is not monotonic") + } + } + if (ID_SHIFT(id) != 0.) + call aaddkd (FITDATA(id,1), ID_SHIFT(id), FITDATA(id,1),ID_NPTS(id)) +end + + +# ID_FITFEATURES -- Compute fit coordinates for features. + +procedure id_fitfeatures (id) + +pointer id # ID pointer +int i + +double id_fitpt() + +begin + if (ID_NFEATURES(id) < 1) + return + + if (ID_CV(id) == NULL) + do i = 1, ID_NFEATURES(id) + FIT(id,i) = id_fitpt (id, PIX(id,i)) + else { + call dcvvector (ID_CV(id), PIX(id,1), FIT(id,1), ID_NFEATURES(id)) + if (ID_SHIFT(id) != 0.) + call aaddkd (FIT(id,1), ID_SHIFT(id), FIT(id,1), ID_NFEATURES(id)) + } +end + + +# ID_FITPT -- Compute fit coordinates from pixel coordinates. + +double procedure id_fitpt (id, pix) + +pointer id # ID pointer +double pix # Pixel coordinate + +double fit + +double smw_c1trand(), shdr_lw(), dcveval() + +begin + if (ID_CV(id) == NULL) { + fit = smw_c1trand (ID_PL(id), pix) + fit = shdr_lw (ID_SH(id), fit) + } else + fit = dcveval (ID_CV(id), pix) + fit = fit + ID_SHIFT(id) + + return (fit) +end + + +# FIT_TO_PIX -- Transform fit coordinate to pixel coordinate. + +define DXMIN .01 + +double procedure fit_to_pix (id, fitcoord) + +pointer id # ID pointer +double fitcoord # Fit coordinate to be transformed +double pixcoord # Pixel coordinate returned + +int i, np1 +double dx + +double smw_c1trand(), id_fitpt() + +begin + np1 = NP1(ID_SH(id)) - 1 + if (FITDATA(id,1) < FITDATA(id,ID_NPTS(id))) { + if ((fitcoord<FITDATA(id,1)) || (fitcoord>FITDATA(id,ID_NPTS(id)))) + return (INDEFD) + + for (i = 1; fitcoord > FITDATA(id,i); i = i + 1) + ; + + if (FITDATA(id,i) == fitcoord) + return (PIXDATA(id,i)) + + pixcoord = smw_c1trand (ID_LP(id), double(i+np1-.5)) + dx = smw_c1trand (ID_LP(id), double(i+np1+.5)) - pixcoord + while (dx > DXMIN) { + dx = dx / 2 + if (id_fitpt (id, pixcoord) < fitcoord) + pixcoord = pixcoord + dx + else + pixcoord = pixcoord - dx + } + } else { + if ((fitcoord<FITDATA(id,ID_NPTS(id))) || (fitcoord>FITDATA(id,1))) + return (INDEFD) + + for (i = 1; fitcoord < FITDATA(id,i); i = i + 1) + ; + + if (FITDATA(id,i) == fitcoord) + return (PIXDATA(id,i)) + + pixcoord = smw_c1trand (ID_LP(id), double(i+np1-.5)) + dx = smw_c1trand (ID_LP(id), double(i+np1+.5)) - pixcoord + while (dx > DXMIN) { + dx = dx / 2 + if (id_fitpt (id, pixcoord) < fitcoord) + pixcoord = pixcoord - dx + else + pixcoord = pixcoord + dx + } + } + + return (pixcoord) +end diff --git a/noao/rv/rvidlines/idfixx.x b/noao/rv/rvidlines/idfixx.x new file mode 100644 index 00000000..ae74fcdf --- /dev/null +++ b/noao/rv/rvidlines/idfixx.x @@ -0,0 +1,27 @@ +include <smw.h> + +# ID_FIXX - Adjust so that pixel indices are increasing. + +procedure id_fixx (sh, x1, x2, y1, y2, i1, i2) + +pointer sh +real x1, x2, y1, y2 +int i1, i2 + +double z, z1, z2, shdr_wl(), shdr_lw() + +begin + z1 = x1 + z2 = x2 + z1 = max (0.5D0, min (double (SN(sh)+.499), shdr_wl(sh, z1))) + z2 = max (0.5D0, min (double (SN(sh)+.499), shdr_wl(sh, z2))) + if (z1 > z2) { + z = y1; y1 = y2; y2 = z + z = z1; z1 = z2; z2 = z + } + + x1 = shdr_lw (sh, z1) + x2 = shdr_lw (sh, z2) + i1 = nint (z1) + i2 = nint (z2) +end diff --git a/noao/rv/rvidlines/idgdata.x b/noao/rv/rvidlines/idgdata.x new file mode 100644 index 00000000..caddce02 --- /dev/null +++ b/noao/rv/rvidlines/idgdata.x @@ -0,0 +1,74 @@ +include <imhdr.h> +include <imio.h> +include <pkg/gtools.h> +include <smw.h> +include "identify.h" + +define SZ_TITLE 320 # Size of long string for title. + +# ID_GDATA -- Get image data. + +procedure id_gdata (id) + +pointer id # ID pointer + +int i, np1 +double hjd +pointer sp, str, im, mw, sh + +double smw_c1trand() +errchk shdr_open, id_vhelio + +begin + call smark (sp) + call salloc (str, SZ_TITLE, TY_CHAR) + + sh = ID_SH(id) + im = IM(sh) + mw = MW(sh) + + # If format is multispec then header info depends on line. + if (SMW_FORMAT(mw) == SMW_ES || SMW_FORMAT(mw) == SMW_MS) + ID_LINE(id,2) = 1 + call shdr_open (im, mw, ID_LINE(id,1), ID_LINE(id,2), + INDEFI, SHDATA, sh) + ID_AP(id,1) = AP(sh) + ID_AP(id,2) = ID_LINE(id,2) + ID_NPTS(id) = SN(sh) + call id_dbsection (id, Memc[ID_IMAGE(id)], ID_AP(id,1), + Memc[ID_SECTION(id)], SZ_FNAME) + call sprintf (Memc[str], SZ_TITLE, "%s %s%s\n%s") + if (ID_TASK(id) == IDENTIFY) + call pargstr ("identify") + else + call pargstr ("rvidlines") + call pargstr (Memc[ID_IMAGE(id)]) + call pargstr (Memc[ID_SECTION(id)]) + call pargstr (TITLE(sh)) + call gt_sets (ID_GT(id), GTTITLE, Memc[str]) + + # Free previous vectors and allocate new vectors. + call mfree (ID_PIXDATA(id), TY_DOUBLE) + + call malloc (ID_PIXDATA(id), ID_NPTS(id), TY_DOUBLE) + + # Set the physical coordinates. + np1 = NP1(sh) - 1 + do i = 1, ID_NPTS(id) + PIXDATA(id,i) = smw_c1trand (ID_LP(id), double(i+np1)) + + # Set the image data + ID_IMDATA(id) = SY(sh) + + # Set the heliocentric correction. + if (ID_TASK(id) == RVIDLINES) { + call id_vhelio (im, ID_ZHELIO(id), hjd, NULL) + ID_ZHELIO(id) = ID_ZHELIO(id) / VLIGHT + } else + ID_ZHELIO(id) = 0 + + ID_NEWGRAPH(id) = YES + ID_NEWCV(id) = YES + + call sfree (sp) +end diff --git a/noao/rv/rvidlines/idgraph.x b/noao/rv/rvidlines/idgraph.x new file mode 100644 index 00000000..c8d92b79 --- /dev/null +++ b/noao/rv/rvidlines/idgraph.x @@ -0,0 +1,168 @@ +include <gset.h> +include <pkg/gtools.h> +include <smw.h> +include "identify.h" + +# ID_GRAPH -- Graph image vector in which features are to be identified. + +procedure id_graph (id, gtype) + +pointer id # ID pointer +int gtype # Graph type + +begin + switch (gtype) { + case 1: + call id_graph1 (id) + case 2: + call id_graph2 (id) + default: + call id_graph1 (id) + } +end + + +procedure id_graph1 (id) + +pointer id # ID pointer + +int i, n +real xmin, xmax, ymin, ymax, dy, xminz, xmaxz, id_zshiftr() +pointer gp, sh, sp, x, y, str + +begin + gp = ID_GP(id) + sh = ID_SH(id) + + call smark (sp) + call salloc (x, SN(sh), TY_REAL) + y = SY(sh) + n = SN(sh) + + call achtdr (FITDATA(id,1), Memr[x], n) + + call gclear (gp) + xmin = min (Memr[x], Memr[x+n-1]) + xmax = max (Memr[x], Memr[x+n-1]) + call alimr (Memr[y], n, ymin, ymax) + dy = ymax - ymin + call gswind (gp, xmin, xmax, ymin - .2 * dy, ymax + .2 * dy) + call gt_swind (gp, ID_GT(id)) + + if (ID_TASK(id) == RVIDLINES && ID_REDSHIFT(id) != 0D0) { + call salloc (str, SZ_LINE, TY_CHAR) + if (ID_ZHELIO(id) == 0D0) + call sprintf (Memc[str], SZ_LINE, + "Vobs = %.5g, Zobs = %.5g\n\n\n") + else + call sprintf (Memc[str], SZ_LINE, + "Vhelio = %.5g, Zhelio = %.5g\n\n\n") + call pargd ((ID_REDSHIFT(id)+ID_ZHELIO(id)) * VLIGHT) + call pargd (ID_REDSHIFT(id)+ID_ZHELIO(id)) + call gt_sets (ID_GT(id), GTSUBTITLE, Memc[str]) + call gseti (gp, G_XDRAWAXES, 1) + call gt_labax (gp, ID_GT(id)) + call gt_sets (ID_GT(id), GTSUBTITLE, "") + + call ggwind (gp, xmin, xmax, ymin, ymax) + xminz = id_zshiftr (id, xmin, 0) + xmaxz = id_zshiftr (id, xmax, 0) + call gswind (gp, xminz, xmaxz, ymin, ymax) + call gseti (gp, G_XDRAWAXES, 2) + call gseti (gp, G_YDRAWAXES, 0) + call glabax (gp, "", "", "") + + call gswind (gp, xmin, xmax, ymin, ymax) + call gctran (gp, xmin, ymin, xmax, ymax, 1, 0) + call gctran (gp, xmax, ymax, xmin, ymin, 0, 1) + } else + call gt_labax (gp, ID_GT(id)) + + call gt_plot (gp, ID_GT(id), Memr[x], Memr[y], n) + + do i = 1, ID_NFEATURES(id) + call id_mark (id, i) + + call sfree (sp) +end + + +# ID_GRAPH2 -- Make review graph for current feature. + +procedure id_graph2 (id) + +pointer id # ID pointer + +int i, j, k, n +real xmin, xmax, ymin, ymax, dy, xminz, xmaxz, id_zshiftr() +pointer gp, sh, sp, x, y, str + +begin + gp = ID_GP(id) + sh = ID_SH(id) + + call smark (sp) + call salloc (x, SN(sh), TY_REAL) + y = SY(sh) + n = SN(sh) + + call achtdr (FITDATA(id,1), Memr[x], n) + + xmin = real (FIT(id,ID_CURRENT(id))) - ID_ZWIDTH(id) / 2. + xmax = real (FIT(id,ID_CURRENT(id))) + ID_ZWIDTH(id) / 2. + + i = 0 + do k = 1, n { + if ((Memr[x+k-1] < xmin) || (Memr[x+k-1] > xmax)) + next + if (i == 0) + i = k + j = k + } + k = j - i + 1 + + call alimr (Memr[y+i-1], k, ymin, ymax) + dy = ymax - ymin + + call gclear (gp) + call gswind (gp, xmin, xmax, ymin - .2 * dy, ymax + .2 * dy) +# if (ID_GT(id) != NULL) { +# call gseti (gp, G_XTRAN, GT_XTRAN(ID_GT(id))) +# call gseti (gp, G_YTRAN, GT_YTRAN(ID_GT(id))) +# } + if (ID_TASK(id) == RVIDLINES && ID_REDSHIFT(id) != 0D0) { + call salloc (str, SZ_LINE, TY_CHAR) + if (ID_ZHELIO(id) == 0D0) + call sprintf (Memc[str], SZ_LINE, + "Vobs = %.5g, Zobs = %.5g\n\n\n") + else + call sprintf (Memc[str], SZ_LINE, + "Vhelio = %.5g, Zhelio = %.5g\n\n\n") + call pargd ((ID_REDSHIFT(id)+ID_ZHELIO(id)) * VLIGHT) + call pargd (ID_REDSHIFT(id)+ID_ZHELIO(id)) + call gt_sets (ID_GT(id), GTSUBTITLE, Memc[str]) + call gseti (gp, G_XDRAWAXES, 1) + call gt_labax (gp, ID_GT(id)) + call gt_sets (ID_GT(id), GTSUBTITLE, "") + + call ggwind (gp, xmin, xmax, ymin, ymax) + xminz = id_zshiftr (id, xmin, 0) + xmaxz = id_zshiftr (id, xmax, 0) + call gswind (gp, xminz, xmaxz, ymin, ymax) + call gseti (gp, G_XDRAWAXES, 2) + call gseti (gp, G_YDRAWAXES, 0) + call glabax (gp, "", "", "") + + call gswind (gp, xmin, xmax, ymin, ymax) + call gctran (gp, xmin, ymin, xmax, ymax, 1, 0) + call gctran (gp, xmax, ymax, xmin, ymin, 0, 1) + } else + call gt_labax (gp, ID_GT(id)) + + call gt_plot (gp, ID_GT(id), Memr[x], Memr[y], n) + + do i = 1, ID_NFEATURES(id) + call id_mark (id, i) + + call sfree (sp) +end diff --git a/noao/rv/rvidlines/ididentify.x b/noao/rv/rvidlines/ididentify.x new file mode 100644 index 00000000..9874865e --- /dev/null +++ b/noao/rv/rvidlines/ididentify.x @@ -0,0 +1,795 @@ +include <error.h> +include <imhdr.h> +include <gset.h> +include <smw.h> +include "identify.h" + +define HELP "noao$onedspec/identify/identify.key" +define RVHELP "noao$rv/rvidlines/rvidlines.key" +define ICFITHELP "noao$lib/scr/idicgfit.key" +define PROMPT "identify options" + +define PAN 1 # Pan graph +define ZOOM 2 # Zoom graph + +# ID_IDENTIFY -- Identify features in an image. +# This is the main interactive loop. + +procedure id_identify (id) + +pointer id # ID pointer + +real wx, wy, wx2, wy2 +int wcs, key +char cmd[SZ_LINE] + +char newimage[SZ_FNAME] +int i, j, last, all, prfeature, nfeatures1, npeaks, ftype, newline[2] +bool answer +double pix, fit, user, shift, pix_shift, z_shift, xg[10] +pointer peaks, label + +bool clgetb() +pointer gopen() +int clgcur(), scan(), nscan(), find_peaks(), errcode(), id_gid(), nowhite() +double id_center(), fit_to_pix(), id_fitpt() +double id_shift(), id_rms(), id_zshiftd(), id_zval() +errchk id_gdata(), id_graph(), id_dbread(), xt_mk1d(), id_log() + +define newim_ 10 +define newkey_ 20 +define beep_ 99 + +begin +newim_ + # Open the image and return if there is an error. + iferr (call id_map (id)) { + call erract (EA_WARN) + return + } + + # Get the image data and return if there is an error. + iferr (call id_gdata (id)) { + call erract (EA_WARN) + return + } + + # Get the database entry for the image if it exists. + iferr { + call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO, YES) + ID_NEWDBENTRY(id) = NO + } then + if ((ID_NFEATURES(id) > 0) || (ID_CV(id) != NULL)) + ID_NEWDBENTRY(id) = YES + + # Set the coordinate information. + iferr (call id_fitdata (id)) + ; + + # Set fitting limits. + call ic_putr (ID_IC(id), "xmin", real (PIXDATA(id,1))) + call ic_putr (ID_IC(id), "xmax", real (PIXDATA(id,ID_NPTS(id)))) + call ic_pstr (ID_IC(id), "help", ICFITHELP) + + # Open graphics. + call clgstr ("graphics", newimage, SZ_FNAME) + ID_GP(id) = gopen (newimage, NEW_FILE, STDGRAPH) + + # Initialize. + ID_GTYPE(id) = PAN + all = 0 + last = ID_CURRENT(id) + newimage[1] = EOS + newline[1] = ID_LINE(id,1) + newline[2] = ID_LINE(id,2) + ID_REFIT(id) = NO + ID_NEWFEATURES(id) = NO + ID_NEWCV(id) = NO + wy = INDEF + key = 'r' + + repeat { + prfeature = YES + if (all != 0) + all = mod (all + 1, 3) + + switch (key) { + case '?': # Print help + if (ID_TASK(id) == IDENTIFY) + call gpagefile (ID_GP(id), HELP, PROMPT) + else + call gpagefile (ID_GP(id), RVHELP, PROMPT) + case ':': # Process colon commands + if (cmd[1] == '/') + call gt_colon (cmd, ID_GP(id), ID_GT(id), ID_NEWGRAPH(id)) + else + call id_colon (id, cmd, newimage, prfeature) + case ' ': # Go to current feature + case '.': # Go to nearest feature + if (ID_NFEATURES(id) == 0) + goto beep_ + call id_nearest (id, double (wx)) + case '-': # Go to previous feature + if (ID_CURRENT(id) == 1) + goto beep_ + ID_CURRENT(id) = ID_CURRENT(id) - 1 + case '+', 'n': # Go to next feature + if (ID_CURRENT(id) == ID_NFEATURES(id)) + goto beep_ + ID_CURRENT(id) = ID_CURRENT(id) + 1 + case 'a': # Set all flag for next key + all = 1 + case 'b': # Mark blended features + i = 1 + fit = wx + xg[i] = fit_to_pix (id, fit) + call printf ("mark other components (exit with 'q'):") + while (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)!=EOF) { + if (key == 'q') + break + i = i + 1 + fit = wx + xg[i] = fit_to_pix (id, fit) + if (i == 10) + break + } + + call printf ("mark two background points: ") + j = clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) + j = clgcur ("cursor", wx2, wy2, wcs, key, cmd, SZ_LINE) + wx = fit_to_pix (id, double (wx)) + wx2 = fit_to_pix (id, double (wx2)) + + switch (ID_FTYPE(id)) { + case EMISSION: + ftype = GEMISSION + case ABSORPTION: + ftype = GABSORPTION + default: + ftype = ID_FTYPE(id) + } + iferr (call id_gcenter (id, xg, i, wx, wy, wx2, wy2, + ID_FWIDTH(id), ftype, YES)) { + call erract (EA_WARN) + prfeature = NO + goto newkey_ + } + + do j = 1, i { + pix = xg[j] + fit = id_fitpt (id, pix) + user = fit + call id_newfeature (id, pix, fit, user, 1.0D0, + ID_FWIDTH(id), ftype, NULL) + USER(id,ID_CURRENT(id)) = INDEFD + call id_match (id, FIT(id,ID_CURRENT(id)), + USER(id,ID_CURRENT(id)), + Memi[ID_LABEL(id)+ID_CURRENT(id)-1], + ID_MATCH(id)) + call id_mark (id, ID_CURRENT(id)) + call printf ("%10.2f %10.8g ") + call pargd (PIX(id,ID_CURRENT(id))) + call pargd (FIT(id,ID_CURRENT(id))) + if (ID_REDSHIFT(id) != 0.) { + call printf ("%10.8g ") + call pargd ( + id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0)) + } + call printf ("(%10.8g %s): ") + call pargd (USER(id,ID_CURRENT(id))) + label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + call gargwrd (cmd, SZ_LINE) + i = nscan() + if (i > 0) { + USER(id,ID_CURRENT(id)) = user + fit = id_zshiftd (id, user, 1) + call id_match (id, fit, USER(id,ID_CURRENT(id)), + Memi[ID_LABEL(id)+ID_CURRENT(id)-1], + ID_MATCH(id)) + } + if (i > 1) { + call reset_scan () + call gargd (user) + call gargstr (cmd, SZ_LINE) + call id_label (cmd, + Memi[ID_LABEL(id)+ID_CURRENT(id)-1]) + } + } + } + case 'c': # Recenter features + if (all != 0) { + for (i = 1; i <= ID_NFEATURES(id); i = i + 1) { + call gseti (ID_GP(id), G_PLTYPE, 0) + call id_mark (id, i) + call gseti (ID_GP(id), G_PLTYPE, 1) + FWIDTH(id,i) = ID_FWIDTH(id) + PIX(id,i) = id_center (id, PIX(id,i), 1, FWIDTH(id,i), + FTYPE(id,i), NO) + if (!IS_INDEFD (PIX(id,i))) { + FIT(id,i) = id_fitpt (id, PIX(id,i)) + call id_mark (id, i) + } else { + call id_delete (id, i) + i = i - 1 + } + } + ID_NEWFEATURES(id) = YES + } else { + if (ID_NFEATURES(id) < 1) + goto beep_ + call id_nearest (id, double (wx)) + pix = PIX(id,ID_CURRENT(id)) + pix = id_center (id, pix, 1, ID_FWIDTH(id), + FTYPE(id,ID_CURRENT(id)), NO) + if (!IS_INDEFD (pix)) { + call gseti (ID_GP(id), G_PLTYPE, 0) + call id_mark (id, ID_CURRENT(id)) + PIX(id,ID_CURRENT(id)) = pix + FWIDTH(id,ID_CURRENT(id)) = ID_FWIDTH(id) + FIT(id,ID_CURRENT(id)) = id_fitpt (id, pix) + call gseti (ID_GP(id), G_PLTYPE, 1) + call id_mark (id, ID_CURRENT(id)) + ID_NEWFEATURES(id) = YES + } else { + call printf ("Centering failed\n") + prfeature = NO + } + } + case 'd': # Delete features + if (all != 0) { + ID_NFEATURES(id) = 0 + ID_CURRENT(id) = 0 + ID_NEWFEATURES(id) = YES + ID_NEWGRAPH(id) = YES + } else { + if (ID_NFEATURES(id) < 1) + goto beep_ + call id_nearest (id, double (wx)) + call gseti (ID_GP(id), G_PLTYPE, 0) + call id_mark (id, ID_CURRENT(id)) + call gseti (ID_GP(id), G_PLTYPE, 1) + call id_delete (id, ID_CURRENT(id)) + ID_CURRENT(id) = min (ID_NFEATURES(id), ID_CURRENT(id)) + last = 0 + } + case 'f': # Fit dispersion function + if (ID_TASK(id) == IDENTIFY) { + call id_dofit (id, YES) + } else { + call id_velocity (id, YES) + prfeature = NO + } + case 'g': # Fit shift + if (ID_TASK(id) == RVIDLINES) + goto beep_ + + call id_doshift (id, YES) + prfeature = NO + case 'i': # Initialize + call dcvfree (ID_CV(id)) + ID_SHIFT(id) = 0. + ID_REDSHIFT(id) = 0. + ID_NEWCV(id) = YES + ID_NFEATURES(id) = 0 + ID_CURRENT(id) = 0 + ID_NEWFEATURES(id) = YES + ID_NEWGRAPH(id) = YES + case 'j': # Go to previous line + newline[1] = ID_LINE(id,1) - ID_NSUM(id,1) + if (newline[1] < 1) { + newline[1] = newline[1] + ID_MAXLINE(id,1) + newline[2] = ID_LINE(id,2) - ID_NSUM(id,2) + if (newline[2] < 1) + newline[2] = newline[2] + ID_MAXLINE(id,2) + } + case 'k': # Go to next line + newline[1] = ID_LINE(id,1) + ID_NSUM(id,1) + if (newline[1] > ID_MAXLINE(id,1)) { + newline[1] = newline[1] - ID_MAXLINE(id,1) + newline[2] = ID_LINE(id,2) + ID_NSUM(id,2) + if (newline[2] > ID_MAXLINE(id,2)) + newline[2] = newline[2] - ID_MAXLINE(id,2) + } + case 'l': # Find features from line list + if (ID_TASK(id) == IDENTIFY) { + if (ID_NFEATURES(id) >= 2) + call id_dofit (id, NO) + if (ID_NEWCV(id) == YES) { + iferr (call id_fitdata(id)) + ; + call id_fitfeatures(id) + ID_NEWCV(id) = NO + } + call id_linelist (id) + if (ID_NEWFEATURES(id) == YES) + ID_REFIT(id) = YES + } else { + call id_velocity (id, NO) + call id_linelist (id) + if (ID_NEWFEATURES(id) == YES) { + call id_velocity (id, NO) + ID_NEWGRAPH(id) = YES + } + } + case 'm': # Mark new feature + fit = wx + pix = fit_to_pix (id, fit) + pix = id_center (id, pix, 1, ID_FWIDTH(id), ID_FTYPE(id), YES) + if (IS_INDEFD (pix)) + goto beep_ + fit = id_fitpt (id, pix) + user = fit + call id_newfeature (id, pix, fit, user, 1.0D0, ID_FWIDTH(id), + ID_FTYPE(id), NULL) + USER(id,ID_CURRENT(id)) = INDEFD + call id_match (id, FIT(id,ID_CURRENT(id)), + USER(id,ID_CURRENT(id)), + Memi[ID_LABEL(id)+ID_CURRENT(id)-1], + ID_MATCH(id)) + call id_mark (id, ID_CURRENT(id)) + call printf ("%10.2f %10.8g ") + call pargd (PIX(id,ID_CURRENT(id))) + call pargd (FIT(id,ID_CURRENT(id))) + if (ID_REDSHIFT(id) != 0.) { + call printf ("%10.8g ") + call pargd ( + id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0)) + } + call printf ("(%10.8g %s): ") + call pargd (USER(id,ID_CURRENT(id))) + label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + call gargwrd (cmd, SZ_LINE) + i = nscan() + if (i > 0) { + USER(id,ID_CURRENT(id)) = user + fit = id_zshiftd (id, user, 1) + call id_match (id, fit, USER(id,ID_CURRENT(id)), + Memi[ID_LABEL(id)+ID_CURRENT(id)-1], + ID_MATCH(id)) + } + if (i > 1) { + call reset_scan () + call gargd (user) + call gargstr (cmd, SZ_LINE) + call id_label (cmd, Memi[ID_LABEL(id)+ID_CURRENT(id)-1]) + } + } + case 'o': # Go to a specified line + call printf ("Line/Column/Band (%d %d): ") + call pargi (ID_LINE(id,1)) + call pargi (ID_LINE(id,2)) + call flush (STDOUT) + if (scan() != EOF) { + call gargi (j) + if (nscan() == 1) { + if (j < 1 || j > ID_MAXLINE(id,1)) + goto beep_ + newline[1] = j + call gargi (j) + if (nscan() == 2) { + if (j < 1 || j > ID_MAXLINE(id,2)) + goto beep_ + newline[2] = j + } + } + } + case 'p': # Switch to pan mode + if (ID_GTYPE(id) != PAN) { + ID_GTYPE(id) = PAN + ID_NEWGRAPH(id) = YES + } + case 'q': # Exit loop + break + case 'r': # Redraw the graph + ID_NEWGRAPH(id) = YES + case 's', 'x': # Shift or correlate features + if (ID_TASK(id) == RVIDLINES) + goto beep_ + + # Get coordinate shift. + switch (key) { + case 's': + call printf ("User coordinate (%10.8g): ") + call pargr (wx) + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + if (nscan() == 1) + shift = wx - user + } else + shift = 0. + case 'x': + if (ID_NFEATURES(id) > 5) + shift = id_shift (id) + else + goto beep_ + } + + ID_NEWFEATURES(id) = YES + ID_NEWCV(id) = YES + ID_NEWGRAPH(id) = YES + prfeature = NO + + if (ID_NFEATURES(id) < 1) { + call printf ("User coordinate shift=%5f\n") + call pargd (shift) + ID_SHIFT(id) = ID_SHIFT(id) + shift + goto newkey_ + } + + # Recenter features. + pix_shift = 0. + z_shift = 0. + nfeatures1 = ID_NFEATURES(id) + + j = 0. + do i = 1, ID_NFEATURES(id) { + pix = fit_to_pix (id, FIT(id,i) + shift) + pix = id_center (id, pix, 1, FWIDTH(id,i), FTYPE(id,i), NO) + if (IS_INDEFD (pix)) { + if (ID_CURRENT(id) == i) + ID_CURRENT(id) = i + 1 + next + } + fit = id_fitpt (id, pix) + + pix_shift = pix_shift + pix - PIX(id,i) + if (FIT(id,i) != 0.) + z_shift = z_shift + id_zval (id, fit, FIT(id,i)) + + j = j + 1 + PIX(id,j) = pix + FIT(id,j) = FIT(id,i) + USER(id,j) = USER(id,i) + WTS(id,j) = WTS(id,i) + FWIDTH(id,j) = FWIDTH(id,i) + FTYPE(id,j) = FTYPE(id,i) + if (ID_CURRENT(id) == i) + ID_CURRENT(id) = j + } + if (j != ID_NFEATURES(id)) { + ID_NFEATURES(id) = j + ID_CURRENT(id) = min (ID_CURRENT(id), ID_NFEATURES(id)) + } + + if (ID_NFEATURES(id) < 1) { + call printf ("User coordinate shift=%5f") + call pargd (shift) + call printf (", No features found during recentering\n") + ID_SHIFT(id) = ID_SHIFT(id) + shift + goto newkey_ + } + + # Adjust shift. + pix = ID_SHIFT(id) + call id_doshift (id, NO) + call id_fitfeatures (id) + + # Print results. + call printf ("Recentered=%d/%d") + call pargi (ID_NFEATURES(id)) + call pargi (nfeatures1) + call printf ( + ", pixel shift=%.2f, user shift=%5f, z=%7.3g, rms=%5g\n") + call pargd (pix_shift / ID_NFEATURES(id)) + call pargd (pix - ID_SHIFT(id)) + call pargd (z_shift / ID_NFEATURES(id)) + call pargd (id_rms(id)) + case 't': # Move the current feature + if (ID_CURRENT(id) < 1) + goto beep_ + pix = fit_to_pix (id, double (wx)) + call gseti (ID_GP(id), G_PLTYPE, 0) + call id_mark (id, ID_CURRENT(id)) + PIX(id,ID_CURRENT(id)) = pix + FIT(id,ID_CURRENT(id)) = id_fitpt (id, pix) + call gseti (ID_GP(id), G_PLTYPE, 1) + call id_mark (id, ID_CURRENT(id)) + ID_NEWFEATURES(id) = YES + case 'u': # Set user coordinate + if (ID_NFEATURES(id) < 1) + goto beep_ + call printf ("%10.2f %10.8g ") + call pargd (PIX(id,ID_CURRENT(id))) + call pargd (FIT(id,ID_CURRENT(id))) + if (ID_REDSHIFT(id) != 0.) { + call printf ("%10.8g ") + call pargd (id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0)) + } + call printf ("(%10.8g %s): ") + call pargd (USER(id,ID_CURRENT(id))) + label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + call gargwrd (cmd, SZ_LINE) + i = nscan() + if (i > 0) { + USER(id,ID_CURRENT(id)) = user + ID_NEWFEATURES(id) = YES + } + if (i > 1) { + call reset_scan () + call gargd (user) + call gargstr (cmd, SZ_LINE) + call id_label (cmd, Memi[ID_LABEL(id)+ID_CURRENT(id)-1]) + } + } + call printf ("Weight (%g): ") + call pargd (WTS(id,ID_CURRENT(id))) + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + if (nscan() > 0) + WTS(id,ID_CURRENT(id)) = user + } + case 'w': # Window graph + call gt_window (ID_GT(id), ID_GP(id), "cursor", ID_NEWGRAPH(id)) + case 'y': # Find peaks + call alimr (IMDATA(id,1), ID_NPTS(id), wx, wy) + call malloc (peaks, ID_NPTS(id), TY_REAL) + if (ID_FTYPE(id) == ABSORPTION || ID_FTYPE(id) == GABSORPTION) + npeaks = find_peaks (IMDATA(id,1), Memr[peaks], + ID_NPTS(id), -1., int (ID_MINSEP(id)), 0, + ID_MAXFEATURES(id), wy, false) + else + npeaks = find_peaks (IMDATA(id,1), Memr[peaks], + ID_NPTS(id), 0., int (ID_MINSEP(id)), 0, + ID_MAXFEATURES(id), wx, false) + for (j = 1; j <= ID_NFEATURES(id); j = j + 1) { + for (i = 1; i <= npeaks; i = i + 1) { + if (!IS_INDEF (Memr[peaks+i-1])) { + pix = Memr[peaks+i-1] + if (abs (pix - PIX(id,j)) < ID_MINSEP(id)) + Memr[peaks+i-1] = INDEF + } + } + } + for (i = 1; i <= npeaks; i = i + 1) { + if (IS_INDEF(Memr[peaks+i-1])) + next + pix = Memr[peaks+i-1] + pix = id_center(id, pix, 1, ID_FWIDTH(id), ID_FTYPE(id), NO) + if (IS_INDEFD (pix)) + next + fit = id_fitpt (id, pix) + user = INDEFD + call id_match (id, fit, user, label, ID_MATCH(id)) + call id_newfeature (id, pix, fit, user, 1.0D0, + ID_FWIDTH(id), ID_FTYPE(id), label) + call id_mark (id, ID_CURRENT(id)) + } + call mfree (peaks, TY_REAL) + case 'z': # Go to zoom mode + if (ID_NFEATURES(id) < 1) + goto beep_ + if (ID_GTYPE(id) == PAN) + ID_NEWGRAPH(id) = YES + ID_GTYPE(id) = ZOOM + call id_nearest (id, double (wx)) + case 'I': + call fatal (0, "Interrupt") + default: +beep_ call printf ("\007") + } + +newkey_ + # Set update flag if anything has changed. + if ((ID_NEWFEATURES(id) == YES) || (ID_NEWCV(id) == YES)) + ID_NEWDBENTRY(id) = YES + + # If a new image exit loop, update database, and start over. + if (newimage[1] != EOS) + break + + # If a new line, save features and set new line. + if (newline[1] != ID_LINE(id,1) || newline[2] != ID_LINE(id,2)) { + call id_saveid (id, ID_LINE(id,1)) + ID_LINE(id,1) = newline[1] + ID_LINE(id,2) = newline[2] + call id_gdata (id) + if (id_gid (id, newline) == EOF) { + iferr { + call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), + NO, NO) + ID_NEWDBENTRY(id) = NO + ID_NEWFEATURES(id) = NO + } then + if ((ID_NFEATURES(id) > 0) || (ID_CV(id) != NULL)) + ID_NEWDBENTRY(id) = YES + } + ID_NEWCV(id) = YES + ID_NEWGRAPH(id) = YES + wy = INDEF + } + + # Refit dispersion function + if (ID_REFIT(id) == YES) { + call id_dofit (id, NO) + ID_REFIT(id) = NO + } + + # If there is a new dispersion solution evaluate the coordinates + if (ID_NEWCV(id) == YES) { + iferr (call id_fitdata (id)) + ; + call id_fitfeatures (id) + ID_NEWCV(id) = NO + } + + # Draw new graph in zoom mode if current feature has changed. + if ((ID_GTYPE(id) == ZOOM) && (last != ID_CURRENT(id))) + ID_NEWGRAPH(id) = YES + + # Draw new graph. + if (ID_NEWGRAPH(id) == YES) { + call id_graph (id, ID_GTYPE(id)) + ID_NEWGRAPH(id) = NO + } + + # Set cursor and print status of current feature (unless canceled). + if (ID_CURRENT(id) > 0) { + if (IS_INDEF (wy)) { + i = max (1, min (ID_NPTS(id), int (PIX(id,ID_CURRENT(id))))) + wy = IMDATA(id,i) + } + + call gscur (ID_GP(id), real (FIT(id,ID_CURRENT(id))), wy) + if (errcode() == OK && prfeature == YES) { + i = ID_CURRENT(id) + pix = PIX(id,i) + fit = FIT(id,i) + user = USER(id,i) + if (ID_TASK(id) == IDENTIFY) { + if (IS_INDEFD(user)) + shift = INDEF + else + shift = fit - user + call printf ("%10.2f %10.8g %10.8g %10.3g %s") + call pargd (pix) + call pargd (fit) + call pargd (user) + call pargd (shift) + label = Memi[ID_LABEL(id)+i-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + } else { + if (IS_INDEFD(user)) + shift = INDEF + else { + shift = id_zval (id, fit, user) - ID_REDSHIFT(id) + if (abs (shift) < 0.01) + shift = shift * VLIGHT + } + call printf ("%10.2f %10.8g %10.8g %10.8g %10.3g %s") + call pargd (pix) + call pargd (fit) + call pargd (id_zshiftd (id, fit, 0)) + call pargd (user) + if (IS_INDEFD(user)) + call pargd (INDEFD) + else + call pargd (shift) + label = Memi[ID_LABEL(id)+i-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + } + } + } + + # Print delayed error message + if (errcode() != OK) + call erract (EA_WARN) + + last = ID_CURRENT(id) + } until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) + call gclose (ID_GP(id)) + + # Write velocity logfile. + if (ID_TASK(id) == RVIDLINES) { + call clgstr ("logfile", cmd, SZ_LINE) + if (nowhite (cmd, cmd, SZ_LINE) == 0) + answer = false + else + answer = true + if (answer && !clgetb ("autowrite")) { + call printf ("Write velocity data to the logfile (yes)? ") + call flush (STDOUT) + if (scan() != EOF) + call gargb (answer) + } + if (answer) { + iferr { + if (ID_NID(id) == 0) + call id_log (id, cmd, NULL) + else { + call id_saveid (id, ID_LINE(id,1)) + for (i = 0; i < ID_NID(id); i = i + 1) { + j = Memi[ID_ID(id)+i] + j =id_gid (id, ID_LINE(j,1)) + call id_dbsection (id, Memc[ID_IMAGE(id)], + ID_AP(id,1), Memc[ID_SECTION(id)], + SZ_FNAME) + call id_log (id, cmd, NULL) + } + } + } then + call erract (EA_WARN) + } + } + + # Warn user that feature data is newer than database entry. + if (ID_NEWDBENTRY(id) == YES) + answer = true + else { + answer = false + for (i = 0; i < ID_NID(id); i = i + 1) { + if (ID_NEWDBENTRY(Memi[ID_ID(id)+i]) == YES) { + answer = true + break + } + } + } + if (answer) { + if (!clgetb ("autowrite")) { + call printf ("Write feature data to the database (yes)? ") + call flush (STDOUT) + if (scan() != EOF) + call gargb (answer) + } + if (answer) { + newline[1] = ID_LINE(id,1) + newline[2] = ID_LINE(id,2) + if (ID_NEWDBENTRY(id) == YES) + call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO) + for (i = 0; i < ID_NID(id); i = i + 1) { + j = Memi[ID_ID(id)+i] + if (ID_NEWDBENTRY(j) == YES && + (ID_LINE(j,1)!=newline[1]||ID_LINE(j,2)!=newline[2])) { + j =id_gid (id, ID_LINE(j,1)) + call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), + NO) + } + } + } + } + + call flush (STDOUT) + + # Free image data. + call mfree (ID_PIXDATA(id), TY_DOUBLE) + call mfree (ID_FITDATA(id), TY_DOUBLE) + call id_free1 (id) + + call smw_close (MW(ID_SH(id))) + call imunmap (IM(ID_SH(id))) + call shdr_close (ID_SH(id)) + + # If a new image was requested with colon command start over. + if (newimage[1] != EOS) { + call strcpy (newimage, Memc[ID_IMAGE(id)], SZ_FNAME) + goto newim_ + } +end diff --git a/noao/rv/rvidlines/idinit.x b/noao/rv/rvidlines/idinit.x new file mode 100644 index 00000000..b17f94a1 --- /dev/null +++ b/noao/rv/rvidlines/idinit.x @@ -0,0 +1,352 @@ +include <gset.h> +include <math/curfit.h> +include "identify.h" + +# ID_INIT -- Allocate identify structure + +procedure id_init (id, taskid) + +int taskid #I Task ID +pointer id #O ID pointer + +begin + call calloc (id, LEN_IDSTRUCT, TY_STRUCT) + + ID_TASK(id) = taskid + ID_NALLOC(id) = 20 + ID_NFEATURES(id) = 0 + ID_CURRENT(id) = 0 + ID_NID(id) = 0 + ID_DT(id) = NULL + + call malloc (ID_IMAGE(id), SZ_FNAME, TY_CHAR) + call malloc (ID_SECTION(id), SZ_FNAME, TY_CHAR) + call malloc (ID_DATABASE(id), SZ_FNAME, TY_CHAR) + call malloc (ID_COORDLIST(id), SZ_FNAME, TY_CHAR) + + call malloc (ID_PIX(id), ID_NALLOC(id), TY_DOUBLE) + call malloc (ID_FIT(id), ID_NALLOC(id), TY_DOUBLE) + call malloc (ID_USER(id), ID_NALLOC(id), TY_DOUBLE) + call malloc (ID_WTS(id), ID_NALLOC(id), TY_DOUBLE) + call malloc (ID_FWIDTHS(id), ID_NALLOC(id), TY_REAL) + call malloc (ID_FTYPES(id), ID_NALLOC(id), TY_INT) + call calloc (ID_LABEL(id), ID_NALLOC(id), TY_POINTER) +end + + +# ID_FREE -- Free identify structure. + +procedure id_free (id) + +pointer id # ID pointer + +int i +pointer ptr + +begin + call mfree (ID_IMAGE(id), TY_CHAR) + call mfree (ID_SECTION(id), TY_CHAR) + call mfree (ID_DATABASE(id), TY_CHAR) + call mfree (ID_COORDLIST(id), TY_CHAR) + call mfree (ID_APS(id), TY_INT) + + ptr = ID_LABEL(id) + do i = 1, ID_NALLOC(id) { + call mfree (Memi[ptr], TY_CHAR) + ptr = ptr + 1 + } + + call mfree (ID_PIX(id), TY_DOUBLE) + call mfree (ID_FIT(id), TY_DOUBLE) + call mfree (ID_USER(id), TY_DOUBLE) + call mfree (ID_WTS(id), TY_DOUBLE) + call mfree (ID_FWIDTHS(id), TY_REAL) + call mfree (ID_FTYPES(id), TY_INT) + call mfree (ID_LABEL(id), TY_POINTER) + + if (ID_DT(id) != NULL) + call dtunmap (ID_DT(id)) + call id_free1 (id) + call id_unmapll (id) + call gt_free (ID_GT(id)) + call dcvfree (ID_CV(id)) + call ic_closed (ID_IC(id)) + + call mfree (id, TY_STRUCT) +end + + +# ID_FREE1 -- Free saved identify structures. + +procedure id_free1 (id) + +pointer id # ID pointer + +int i, j +pointer id1, ptr + +begin + for (i = 0; i < ID_NID(id); i = i + 1) { + id1 = Memi[ID_ID(id)+i] + + ptr = ID_LABEL(id1) + do j = 1, ID_NALLOC(id1) { + call mfree (Memi[ptr], TY_CHAR) + ptr = ptr + 1 + } + + call mfree (ID_PIX(id1), TY_DOUBLE) + call mfree (ID_FIT(id1), TY_DOUBLE) + call mfree (ID_USER(id1), TY_DOUBLE) + call mfree (ID_WTS(id1), TY_DOUBLE) + call mfree (ID_FWIDTHS(id1), TY_REAL) + call mfree (ID_FTYPES(id1), TY_INT) + call mfree (ID_LABEL(id1), TY_POINTER) + if (ID_CV(id1) != NULL) + call dcvfree (ID_CV(id1)) + if (ID_IC(id1) != NULL) + call ic_closed (ID_IC(id1)) + call shdr_close (ID_SH(id1)) + call mfree (id1, TY_STRUCT) + } + call mfree (ID_ID(id), TY_POINTER) + + ID_NID(id) = 0 +end + + +# ID_SAVEID -- Save identify information by line. + +procedure id_saveid (id, line) + +pointer id # IDENTIFY structure +int line[2] # Save as line + +int i +pointer id1 + +begin + # Check if already saved. If not saved allocate memory. + for (i = 1; i <= ID_NID(id); i = i + 1) { + id1 = Memi[ID_ID(id)+i-1] + if (ID_LINE(id1,1) == line[1] && ID_LINE(id1,2) == line[2]) + break + } + + call id_sid (id, i) +end + + + +# ID_SID -- Save identify information by index. + +procedure id_sid (id, n) + +pointer id # IDENTIFY structure +int n # Save as index + +int i, j, dcvstati(), strlen() +pointer sp, id1, coeffs, ptr1, ptr2 + +begin + if (n > ID_NID(id)) { + if (n == 1) + call malloc (ID_ID(id), 1, TY_POINTER) + else + call realloc (ID_ID(id), n, TY_POINTER) + call calloc (id1, LEN_IDSTRUCT, TY_STRUCT) + Memi[ID_ID(id)+n-1] = id1 + ID_NID(id) = n + } else + id1 = Memi[ID_ID(id)+n-1] + + # Allocate or reallocate memory for features and copy them. + if (ID_NFEATURES(id) > 0) { + if (ID_NALLOC(id1) == 0) { + call malloc (ID_PIX(id1), ID_NFEATURES(id), TY_DOUBLE) + call malloc (ID_FIT(id1), ID_NFEATURES(id), TY_DOUBLE) + call malloc (ID_USER(id1), ID_NFEATURES(id), TY_DOUBLE) + call malloc (ID_WTS(id1), ID_NFEATURES(id), TY_DOUBLE) + call malloc (ID_FWIDTHS(id1), ID_NFEATURES(id), TY_REAL) + call malloc (ID_FTYPES(id1), ID_NFEATURES(id), TY_INT) + call calloc (ID_LABEL(id1), ID_NFEATURES(id), TY_POINTER) + } else if (ID_NALLOC(id1) != ID_NFEATURES(id)) { + call realloc (ID_PIX(id1), ID_NFEATURES(id), TY_DOUBLE) + call realloc (ID_FIT(id1), ID_NFEATURES(id), TY_DOUBLE) + call realloc (ID_USER(id1), ID_NFEATURES(id), TY_DOUBLE) + call realloc (ID_WTS(id1), ID_NFEATURES(id), TY_DOUBLE) + call realloc (ID_FWIDTHS(id1), ID_NFEATURES(id), TY_REAL) + call realloc (ID_FTYPES(id1), ID_NFEATURES(id), TY_INT) + call realloc (ID_LABEL(id1), ID_NFEATURES(id), TY_POINTER) + + j = ID_NALLOC(id1) + i = ID_NFEATURES(id) - j + if (i > 0) + call aclri (Memi[ID_LABEL(id1)+j], i) + } + call amovd (PIX(id,1), PIX(id1,1), ID_NFEATURES(id)) + call amovd (FIT(id,1), FIT(id1,1), ID_NFEATURES(id)) + call amovd (USER(id,1), USER(id1,1), ID_NFEATURES(id)) + call amovd (WTS(id,1), WTS(id1,1), ID_NFEATURES(id)) + call amovr (FWIDTH(id,1), FWIDTH(id1,1), ID_NFEATURES(id)) + call amovi (FTYPE(id,1), FTYPE(id1,1), ID_NFEATURES(id)) + + ptr1 = ID_LABEL(id) + ptr2 = ID_LABEL(id1) + do i = 1, ID_NFEATURES(id) { + call mfree (Memi[ptr2], TY_CHAR) + if (Memi[ptr1] != NULL) { + j = strlen (Memc[Memi[ptr1]]) + call malloc (Memi[ptr2], j, TY_CHAR) + call strcpy (Memc[Memi[ptr1]], Memc[Memi[ptr2]], j) + } + ptr1 = ptr1 + 1 + ptr2 = ptr2 + 1 + } + + ID_NALLOC(id1) = ID_NFEATURES(id) + } + + # Use a SAVE and RESTORE to copy the CURFIT data. + if (ID_CV(id1) != NULL) + call dcvfree (ID_CV(id1)) + if (ID_CV(id) != NULL) { + call smark (sp) + i = dcvstati (ID_CV(id), CVNSAVE) + call salloc (coeffs, i, TY_DOUBLE) + call dcvsave (ID_CV(id), Memd[coeffs]) + call dcvrestore (ID_CV(id1), Memd[coeffs]) + call sfree (sp) + + if (ID_IC(id1) == NULL) + call ic_open (ID_IC(id1)) + call ic_copy (ID_IC(id), ID_IC(id1)) + } + + #ID_LINE(id1,1) = line[1] + #ID_LINE(id1,2) = line[2] + ID_LINE(id1,1) = ID_LINE(id,1) + ID_LINE(id1,2) = ID_LINE(id,2) + ID_AP(id1,1) = ID_AP(id,1) + ID_AP(id1,2) = ID_AP(id,2) + ID_NFEATURES(id1) = ID_NFEATURES(id) + ID_SHIFT(id1) = ID_SHIFT(id) + ID_REDSHIFT(id1) = ID_REDSHIFT(id) + ID_RMSRED(id1) = ID_RMSRED(id) + ID_ZHELIO(id1) = ID_ZHELIO(id) + ID_CURRENT(id1) = ID_CURRENT(id) + + ID_NEWFEATURES(id1) = ID_NEWFEATURES(id) + ID_NEWCV(id1) = ID_NEWCV(id) + ID_NEWDBENTRY(id1) = ID_NEWDBENTRY(id) +end + + +# ID_GID -- Get saved identify information for specified line. + +int procedure id_gid (id, line) + +pointer id # IDENTIFY structure +int line[2] # Line number to get + +int i, id_getid() + +begin + # Check if saved. + for (i = 1; i <= ID_NID(id); i = i + 1) { + if (ID_LINE(Memi[ID_ID(id)+i-1],1) == line[1] && + ID_LINE(Memi[ID_ID(id)+i-1],2) == line[2]) + break + } + + return (id_getid (id, i)) +end + + +# ID_GETID -- Get saved identify information for specified index. + +int procedure id_getid (id, n) + +pointer id # IDENTIFY structure +int n # Index of saved features to be returned + +int i, j, dcvstati(), strlen() +pointer sp, id1, coeffs, ptr1, ptr2 + +begin + if (n < 1 || n > ID_NID(id)) + return (EOF) + + id1 = Memi[ID_ID(id)+n-1] + + # Reallocate memory for features and copy them. + if (ID_NFEATURES(id1) > 0) { + if (ID_NALLOC(id1) != ID_NALLOC(id)) { + call realloc (ID_PIX(id), ID_NALLOC(id1), TY_DOUBLE) + call realloc (ID_FIT(id), ID_NALLOC(id1), TY_DOUBLE) + call realloc (ID_USER(id), ID_NALLOC(id1), TY_DOUBLE) + call realloc (ID_WTS(id), ID_NALLOC(id1), TY_DOUBLE) + call realloc (ID_FWIDTHS(id), ID_NALLOC(id1), TY_REAL) + call realloc (ID_FTYPES(id), ID_NALLOC(id1), TY_INT) + call realloc (ID_LABEL(id), ID_NALLOC(id1), TY_POINTER) + + j = ID_NALLOC(id) + i = ID_NALLOC(id1) - j + if (i > 0) + call aclri (Memi[ID_LABEL(id)+j], i) + } + call amovd (PIX(id1,1), PIX(id,1), ID_NFEATURES(id1)) + call amovd (FIT(id1,1), FIT(id,1), ID_NFEATURES(id1)) + call amovd (USER(id1,1), USER(id,1), ID_NFEATURES(id1)) + call amovd (WTS(id1,1), WTS(id,1), ID_NFEATURES(id1)) + call amovr (FWIDTH(id1,1), FWIDTH(id,1), ID_NFEATURES(id1)) + call amovi (FTYPE(id1,1), FTYPE(id,1), ID_NFEATURES(id1)) + + ptr1 = ID_LABEL(id1) + ptr2 = ID_LABEL(id) + do i = 1, ID_NFEATURES(id1) { + call mfree (Memi[ptr2], TY_CHAR) + if (Memi[ptr1] != NULL) { + j = strlen (Memc[Memi[ptr1]]) + call malloc (Memi[ptr2], j, TY_CHAR) + call strcpy (Memc[Memi[ptr1]], Memc[Memi[ptr2]], j) + } + ptr1 = ptr1 + 1 + ptr2 = ptr2 + 1 + } + + ID_NALLOC(id) = ID_NALLOC(id1) + ID_NFEATURES(id) = ID_NFEATURES(id1) + ID_NEWFEATURES(id) = ID_NEWFEATURES(id1) + ID_CURRENT(id) = ID_CURRENT(id1) + ID_NEWDBENTRY(id) = ID_NEWDBENTRY(id1) + } + + # Use a SAVE and RESTORE to copy the CURFIT data. + if (ID_CV(id1) != NULL) { + if (ID_CV(id) != NULL) + call dcvfree (ID_CV(id)) + call smark (sp) + i = dcvstati (ID_CV(id1), CVNSAVE) + call salloc (coeffs, i, TY_DOUBLE) + call dcvsave (ID_CV(id1), Memd[coeffs]) + call dcvrestore (ID_CV(id), Memd[coeffs]) + call sfree (sp) + + call ic_copy (ID_IC(id1), ID_IC(id)) + + ID_SHIFT(id) = ID_SHIFT(id1) + ID_REDSHIFT(id) = ID_REDSHIFT(id1) + ID_RMSRED(id) = ID_RMSRED(id1) + ID_ZHELIO(id) = ID_ZHELIO(id1) + ID_NEWCV(id) = ID_NEWCV(id1) + ID_NEWDBENTRY(id) = ID_NEWDBENTRY(id1) + } + + ID_LINE(id,1) = ID_LINE(id1,1) + ID_LINE(id,2) = ID_LINE(id1,2) + ID_AP(id,1) = ID_AP(id1,1) + ID_AP(id,2) = ID_AP(id1,2) + + return (OK) +end diff --git a/noao/rv/rvidlines/idlabel.x b/noao/rv/rvidlines/idlabel.x new file mode 100644 index 00000000..cb5fa439 --- /dev/null +++ b/noao/rv/rvidlines/idlabel.x @@ -0,0 +1,30 @@ +define SKIP ($1==' '||$1=='\t'||$1=='"'||$1=='\'') + +# ID_LABEL -- Set label + +procedure id_label (str, label) + +char str[ARB] # String to be set +pointer label # Label pointer to be set + +int i, j, strlen() +pointer cp + +begin + call mfree (label, TY_CHAR) + + for (i=1; str[i]!=EOS && SKIP(str[i]); i=i+1) + ; + for (j=strlen(str); j>=i && SKIP(str[j]); j=j-1) + ; + + if (i <= j) { + call malloc (label, j-i+1, TY_CHAR) + cp = label + for (; i<=j; i=i+1) { + Memc[cp] = str[i] + cp = cp + 1 + } + Memc[cp] = EOS + } +end diff --git a/noao/rv/rvidlines/idlinelist.x b/noao/rv/rvidlines/idlinelist.x new file mode 100644 index 00000000..d427aad7 --- /dev/null +++ b/noao/rv/rvidlines/idlinelist.x @@ -0,0 +1,250 @@ +include <error.h> +include <mach.h> +include "identify.h" + +# ID_MAPLL -- Read the line list into memory. + +procedure id_mapll (id) + +pointer id # Identify structure + +int fd, nalloc, nlines +pointer ll1, ll2, str + +int open(), fscan(), nscan(), nowhite() +double value + +errchk open, fscan, malloc, realloc + +begin + ID_LL(id) = NULL + + if (nowhite (Memc[ID_COORDLIST(id)], + Memc[ID_COORDLIST(id)], SZ_FNAME) == 0) + return + iferr (fd = open (Memc[ID_COORDLIST(id)], READ_ONLY, TEXT_FILE)) { + call erract (EA_WARN) + return + } + + call malloc (str, SZ_LINE, TY_CHAR) + nalloc = 0 + nlines = 0 + while (fscan (fd) != EOF) { + call gargd (value) + if (nscan() != 1) + next + + if (nalloc == 0) { + nalloc = 100 + call malloc (ll1, nalloc, TY_DOUBLE) + call calloc (ll2, nalloc, TY_POINTER) + } else if (nlines == nalloc) { + nalloc = nalloc + 100 + call realloc (ll1, nalloc, TY_DOUBLE) + call realloc (ll2, nalloc, TY_POINTER) + call aclri (Memi[ll2+nalloc-100], 100) + } + + Memd[ll1+nlines] = value + call gargstr (Memc[str], SZ_LINE) + call id_label (Memc[str], Memi[ll2+nlines]) + + nlines = nlines + 1 + } + call mfree (str, TY_CHAR) + call close (fd) + + if (nlines == 0) + return + + call realloc (ll1, nlines + 1, TY_DOUBLE) + call realloc (ll2, nlines + 1, TY_POINTER) + Memd[ll1+nlines] = INDEFD + call malloc (ID_LL(id), 2, TY_POINTER) + Memi[ID_LL(id)] = ll1 + Memi[ID_LL(id)+1] = ll2 +end + + +# ID_UNMAPLL -- Unmap the linelist. + +procedure id_unmapll (id) + +pointer id # Identify structure + +pointer ll1, ll2 + +begin + if (ID_LL(id) == NULL) + return + + ll1 = Memi[ID_LL(id)] + ll2 = Memi[ID_LL(id)+1] + while (!IS_INDEFD(Memd[ll1])) { + call mfree (Memi[ll2], TY_CHAR) + ll1 = ll1 + 1 + ll2 = ll2 + 1 + } + + call mfree (Memi[ID_LL(id)], TY_DOUBLE) + call mfree (Memi[ID_LL(id)+1], TY_POINTER) + call mfree (ID_LL(id), TY_POINTER) +end + + + +# ID_MATCH -- Match current feature against a line list. +# +# This is extremely inefficient. It can be greatly improved. + +procedure id_match (id, in, out, label, diff) + +pointer id # Identify structure +double in # Coordinate to be matched +double out # Matched coordinate +pointer label # Pointer to label +real diff # Maximum difference + +double zin, delta, deltamin, id_zshiftd() +pointer ll1, ll2, tmp +int strlen() + +begin + if (ID_LL(id) == NULL) { + out = id_zshiftd (id, in, 0) + label = NULL + return + } + + zin = id_zshiftd (id, in, 0) + deltamin = MAX_REAL + + ll1 = Memi[ID_LL(id)] + ll2 = Memi[ID_LL(id)+1] + while (!IS_INDEFD(Memd[ll1])) { + delta = abs (zin - Memd[ll1]) + if (delta < deltamin) { + deltamin = delta + if (deltamin <= diff) { + out = Memd[ll1] + label = Memi[ll2] + } + } + ll1 = ll1 + 1 + ll2 = ll2 + 1 + } + + if (label != NULL) { + call malloc (tmp, strlen (Memc[label]), TY_CHAR) + call strcpy (Memc[label], Memc[tmp], ARB) + label = tmp + } +end + +# ID_LINELIST -- Add features from a line list. + +procedure id_linelist (id) + +pointer id # Identify structure + +int i, nfound, nextpix, lastpix, cursave +double pix, fit, fit1, fit2, user, peak, minval, diff, diff1 +pointer sp, pixes, fits, users, labels, ll1, ll2, label + +double id_center(), fit_to_pix(), id_fitpt(), id_peak(), id_zshiftd() + +begin + if (ID_LL(id) == NULL) + return + + call smark (sp) + call salloc (pixes, ID_MAXFEATURES(id), TY_DOUBLE) + call salloc (fits, ID_MAXFEATURES(id), TY_DOUBLE) + call salloc (users, ID_MAXFEATURES(id), TY_DOUBLE) + call salloc (labels, ID_MAXFEATURES(id), TY_POINTER) + + nfound = 0 + lastpix = 0 + minval = MAX_REAL + + fit1 = min (FITDATA(id,1), FITDATA(id,ID_NPTS(id))) + fit2 = max (FITDATA(id,1), FITDATA(id,ID_NPTS(id))) + ll1 = Memi[ID_LL(id)] + ll2 = Memi[ID_LL(id)+1] + while (!IS_INDEFD(Memd[ll1])) { + user = id_zshiftd (id, Memd[ll1], 1) + label = Memi[ll2] + ll1 = ll1 + 1 + ll2 = ll2 + 1 + if (user < fit1) + next + if (user > fit2) + break + + pix = id_center (id, fit_to_pix (id, user), 1, ID_FWIDTH(id), + ID_FTYPE(id), NO) + if (!IS_INDEFD(pix)) { + fit = id_fitpt (id, pix) + diff = abs (fit - user) + if (diff > ID_MATCH(id)) + next + if (lastpix > 0) { + if (abs (pix - Memd[pixes+lastpix-1]) < 0.01) { + diff1 = abs (Memd[fits+lastpix-1]-Memd[users+lastpix-1]) + if (diff < diff1) { + Memd[pixes+lastpix-1] = pix + Memd[fits+lastpix-1] = fit + Memd[users+lastpix-1] = id_zshiftd (id, user, 0) + Memi[labels+lastpix-1] = label + } + next + } + } + + peak = abs (id_peak (id, pix)) + if (nfound < ID_MAXFEATURES(id)) { + nfound = nfound + 1 + if (peak < minval) { + nextpix = nfound + minval = peak + } + Memd[pixes+nfound-1] = pix + Memd[fits+nfound-1] = id_fitpt (id, pix) + Memd[users+nfound-1] = id_zshiftd (id, user, 0) + Memi[labels+nfound-1] = label + lastpix = nfound + } else if (peak > minval) { + Memd[pixes+nextpix-1] = pix + Memd[fits+nextpix-1] = id_fitpt (id, pix) + Memd[users+nextpix-1] = id_zshiftd (id, user, 0) + Memi[labels+nextpix-1] = label + lastpix = nextpix + + minval = MAX_REAL + do i = 1, nfound { + pix = Memd[pixes+i-1] + peak = abs (id_peak (id, pix)) + if (peak < minval) { + nextpix = i + minval = peak + } + } + } + } + } + + do i = 1, nfound { + pix = Memd[pixes+i-1] + fit = Memd[fits+i-1] + user = Memd[users+i-1] + label = Memi[labels+i-1] + call id_newfeature (id, pix, fit, user, 1.0D0, ID_FWIDTH(id), + ID_FTYPE(id), label) + if (i == 1) + cursave = ID_CURRENT(id) + } + ID_CURRENT(id) = cursave + + call sfree (sp) +end diff --git a/noao/rv/rvidlines/idlog.x b/noao/rv/rvidlines/idlog.x new file mode 100644 index 00000000..58645fc1 --- /dev/null +++ b/noao/rv/rvidlines/idlog.x @@ -0,0 +1,190 @@ +include <smw.h> +include <time.h> +include "identify.h" + + +# ID_LOG -- Write log + +procedure id_log (id, file, logfd) + +pointer id # ID pointer +char file[ARB] # Log file +int logfd # Log fd + +char str[SZ_TIME] +int i, fd, nrms +double z, zrms, zerr, zhelio, resid, rms, v, verr, vhelio, hjd + +int open() +double id_zshiftd(), id_zval() +long clktime() +errchk open, id_velocity, id_vhelio + +begin + if (ID_NFEATURES(id) == 0) + return + + if (logfd == NULL) + fd = open (file, APPEND, TEXT_FILE) + else + fd = logfd + + if (ID_TASK(id) == IDENTIFY) { + + call cnvtime (clktime (0), str, SZ_TIME) + call fprintf (fd, "\n%s\n") + call pargstr (str) + call fprintf (fd, "Features identified in image %s%s: %s\n") + call pargstr (Memc[ID_IMAGE(id)]) + call pargstr (Memc[ID_SECTION(id)]) + call pargstr (TITLE(ID_SH(id))) + + call fprintf (fd, " %8s %10s %10s %10s %5s %s\n") + call pargstr ("Pixel") + call pargstr ("Fit") + call pargstr ("User") + call pargstr ("Residual") + call pargstr ("Fwidth") + call pargstr ("Wt") + call pargstr ("Label") + + rms = 0. + nrms = 0 + do i = 1, ID_NFEATURES(id) { + call fprintf (fd, + "%2d %8.2f %10.8g %10.8g %10.4g %4f %s\n") + call pargi (i) + call pargd (PIX(id,i)) + call pargd (FIT(id,i)) + call pargd (USER(id,i)) + if (IS_INDEFD (USER(id,i))) + call pargd (USER(id,i)) + else { + resid = FIT(id,i) - USER(id,i) + call pargd (resid) + if (WTS(id,i) > 0.) { + rms = rms + resid ** 2 + nrms = nrms + 1 + } + } + call pargd (WTS(id,i)) + if (Memi[ID_LABEL(id)+i-1] != NULL) + call pargstr (Memc[Memi[ID_LABEL(id)+i-1]]) + else + call pargstr ("") + } + + if (nrms > 1) { + call fprintf (fd, "RMS = %0.6g\n") + call pargd (sqrt (rms / nrms)) + } + + } else { + call id_velocity (id, NO) + z = ID_REDSHIFT(id) + zrms = ID_RMSRED(id) + zhelio = ID_ZHELIO(id) + v = ID_REDSHIFT(id) * VLIGHT + call id_vhelio (IM(ID_SH(id)), vhelio, hjd, fd) + + call cnvtime (clktime (0), str, SZ_TIME) + call fprintf (fd, "\n%s\n") + call pargstr (str) + call fprintf (fd, "Features identified in image %s%s: %s\n") + call pargstr (Memc[ID_IMAGE(id)]) + call pargstr (Memc[ID_SECTION(id)]) + call pargstr (TITLE(ID_SH(id))) + + call fprintf (fd, "%10s %10s %10s %10s %10s %5s %s\n") + call pargstr ("Measured") + call pargstr ("User") + call pargstr ("Residual") + call pargstr ("Velocity") + call pargstr ("Residual") + call pargstr ("Wt") + call pargstr ("Label") + + rms = 0. + nrms = 0 + do i = 1, ID_NFEATURES(id) { + call fprintf (fd, + "%10.8g %10.8g %10.4g %10.8g %10.4g %4f %s\n") + call pargd (id_zshiftd (id, FIT(id,i), 0)) + if (IS_INDEFD (USER(id,i))) { + call pargd (INDEFD) + call pargd (INDEFD) + call pargd (INDEFD) + call pargd (INDEFD) + } else { + call pargd (USER(id,i)) + resid = id_zshiftd (id, FIT(id,i), 0) - USER(id,i) + call pargd (resid) + if (WTS(id,i) > 0.) { + rms = rms + resid ** 2 + nrms = nrms + 1 + } + verr = id_zval (id, FIT(id,i), USER(id,i)) * VLIGHT + call pargd (verr + vhelio) + call pargd (verr - v) + } + call pargd (WTS(id,i)) + if (Memi[ID_LABEL(id)+i-1] != NULL) + call pargstr (Memc[Memi[ID_LABEL(id)+i-1]]) + else + call pargstr ("") + } + + if (nrms > 1) { + call fprintf (fd, "Wavelength RMS = %0.6g\n") + call pargd (sqrt (rms / nrms)) + call fprintf (fd, "Velocity RMS = %8.5g\n") + call pargd (zrms * VLIGHT) + } + + zerr = zrms + if (nrms > 1) + zerr = zerr / sqrt (nrms - 1.) + v = z * VLIGHT + verr = zerr * VLIGHT + + call fprintf (fd, "\n") + call fprintf (fd, "%s %3d : Zobs = %10.5g, ") + call pargstr (Memc[ID_IMAGE(id)]) + call pargi (ID_AP(id,1)) + call pargd (z) + call fprintf (fd, "Mean err = %10.5g, Lines = %3d\n") + call pargd (zerr) + call pargi (nrms) + call fprintf (fd, "%s %3d : Vobs = %8.5g km/s, ") + call pargstr (Memc[ID_IMAGE(id)]) + call pargi (ID_AP(id,1)) + call pargd (v) + call fprintf (fd, "Mean err = %8.5g km/s, Lines = %3d\n") + call pargd (verr) + call pargi (nrms) + if (zhelio != 0D0) { + call fprintf (fd, "%s %3d : Zhelio = %10.5g, ") + call pargstr (Memc[ID_IMAGE(id)]) + call pargi (ID_AP(id,1)) + call pargd (z + zhelio) + call fprintf (fd, "Mean err = %8.5g km/s, Lines = %3d\n") + call pargd (zerr) + call pargi (nrms) + call fprintf (fd, "%s %3d : Vhelio = %8.5g km/s, ") + call pargstr (Memc[ID_IMAGE(id)]) + call pargi (ID_AP(id,1)) + call pargd (v + vhelio) + call fprintf (fd, "Mean err = %8.5g km/s, Lines = %3d\n") + call pargd (verr) + call pargi (nrms) + call fprintf (fd, "%s %3d : HJD = %g\n") + call pargstr (Memc[ID_IMAGE(id)]) + call pargi (ID_AP(id,1)) + call pargd (hjd) + } + call fprintf (fd, "\n") + } + + if (logfd == NULL) + call close (fd) +end diff --git a/noao/rv/rvidlines/idmap.x b/noao/rv/rvidlines/idmap.x new file mode 100644 index 00000000..5af57217 --- /dev/null +++ b/noao/rv/rvidlines/idmap.x @@ -0,0 +1,379 @@ +include <ctype.h> +include <imhdr.h> +include <pkg/gtools.h> +include <smw.h> +include "identify.h" + +# Sepcial section words. +define SPECIAL "|first|middle|x|y|z|last|column|line|band|" +define FIRST 1 +define MIDDLE 2 +define X 3 +define Y 4 +define Z 5 +define LAST 6 +define COLUMN 7 +define LINE 8 +define BAND 9 + +# ID_MAP -- Map an image for IDENTIFY/REIDENTIFY +# The image must 1, 2, or 3 dimensional. An image section may be given with +# the image name or with the CL parameter "section". The CL parameter can +# have one of the following formats: +# 1. An IMIO image section +# 2. [line|column|x|y|z] [#|middle|last] [#|middle|last] +# 3. [#|middle|last] [#|middle|last] [line|column|x|y|z] +# where # is a line or column number. The strings may be abbreviated. +# The task returns and error if it cannot map the image or determine +# the 1D line or column desired. + +procedure id_map (id) + +pointer id # IDENTIFY data structure pointer + +int i, j, k, l, a, b, c, x1[3], x2[3], xs[3] +pointer sp, wrd1, wrd2, wrd3, im + +int imaccess(), strdic(), ctoi(), nscan() +pointer immap() +errchk immap, id_maphdr + +begin + # Separate the image name and image section and map the full image. + call imgsection (Memc[ID_IMAGE(id)], Memc[ID_SECTION(id)], SZ_FNAME) + call imgimage (Memc[ID_IMAGE(id)], Memc[ID_IMAGE(id)], SZ_FNAME) + call id_noextn (Memc[ID_IMAGE(id)]) + im = immap (Memc[ID_IMAGE(id)], READ_ONLY, 0) + + # If no image section is found use the "section" parameter. + if (Memc[ID_SECTION(id)] == EOS && IM_NDIM(im) > 1) { + call clgstr ("section", Memc[ID_SECTION(id)], SZ_FNAME) + call xt_stripwhite (Memc[ID_SECTION(id)]) + + # If not an image section construct one. + if (Memc[ID_SECTION(id)] != '[') { + call smark (sp) + call salloc (wrd1, SZ_FNAME, TY_CHAR) + call salloc (wrd2, SZ_FNAME, TY_CHAR) + call salloc (wrd3, SZ_FNAME, TY_CHAR) + + call sscan (Memc[ID_SECTION(id)]) + + # Parse axis and elements. + call gargwrd (Memc[wrd1], SZ_FNAME) + call gargwrd (Memc[wrd2], SZ_FNAME) + call gargwrd (Memc[wrd3], SZ_FNAME) + switch (nscan()) { + case 0: + a = X + b = MIDDLE + c = MIDDLE + case 1: + a = strdic (Memc[wrd1], Memc[wrd1], SZ_FNAME, SPECIAL) + b = MIDDLE + c = MIDDLE + case 2: + a = strdic (Memc[wrd1], Memc[wrd1], SZ_FNAME, SPECIAL) + if (a >= X) + b = strdic (Memc[wrd2], Memc[wrd2], SZ_FNAME, SPECIAL) + else { + b = a + a = strdic (Memc[wrd2], Memc[wrd2], SZ_FNAME, SPECIAL) + call strcpy (Memc[wrd1], Memc[wrd2], SZ_FNAME) + } + c = MIDDLE + call strcpy (Memc[wrd2], Memc[wrd3], SZ_FNAME) + case 3: + a = strdic (Memc[wrd1], Memc[wrd1], SZ_FNAME, SPECIAL) + if (a >= X) { + b = strdic (Memc[wrd2], Memc[wrd2], SZ_FNAME, SPECIAL) + c = strdic (Memc[wrd3], Memc[wrd3], SZ_FNAME, SPECIAL) + } else { + b = a + a = strdic (Memc[wrd2], Memc[wrd2], SZ_FNAME, SPECIAL) + if (a >= X) { + c = strdic (Memc[wrd3], Memc[wrd3],SZ_FNAME,SPECIAL) + call strcpy (Memc[wrd1], Memc[wrd2], SZ_FNAME) + } else { + c = b + b = a + a = strdic (Memc[wrd3], Memc[wrd3],SZ_FNAME,SPECIAL) + call strcpy (Memc[wrd2], Memc[wrd3], SZ_FNAME) + call strcpy (Memc[wrd1], Memc[wrd2], SZ_FNAME) + } + } + } + + switch (a) { + case X, LINE: + i = 1 + j = 2 + k = 3 + case Y, COLUMN: + i = 2 + j = 1 + k = 3 + case Z, BAND: + i = 3 + j = 1 + k = 2 + default: + call imunmap (im) + call error (1, + "Error in section specification or non-unique abbreviation") + } + + switch (b) { + case FIRST: + ID_LINE(id,1) = 1 + case MIDDLE: + ID_LINE(id,1) = (1 + IM_LEN(im,j)) / 2 + case LAST: + ID_LINE(id,1) = IM_LEN(im,j) + default: + l = 1 + if (ctoi (Memc[wrd2], l, ID_LINE(id,1)) == 0) + call error (1, "Error in section specification") + } + + switch (c) { + case FIRST: + ID_LINE(id,2) = 1 + case MIDDLE: + ID_LINE(id,2) = (1 + IM_LEN(im,k)) / 2 + case LAST: + ID_LINE(id,2) = IM_LEN(im,k) + default: + l = 1 + if (ctoi (Memc[wrd3], l, ID_LINE(id,2)) == 0) + call error (1, "Error in section specification") + } + + # Format section. + switch (IM_NDIM(im)) { + case 2: + switch (i) { + case 1: + call sprintf (Memc[ID_SECTION(id)], SZ_FNAME, "[*,%d]") + case 2: + call sprintf (Memc[ID_SECTION(id)], SZ_FNAME, "[%d,*]") + default: + call error (1, "Error in section specification") + } + call pargi (ID_LINE(id,1)) + case 3: + switch (i) { + case 1: + call sprintf (Memc[ID_SECTION(id)],SZ_FNAME,"[*,%d,%d]") + case 2: + call sprintf (Memc[ID_SECTION(id)],SZ_FNAME,"[%d,*,%d]") + case 3: + call sprintf (Memc[ID_SECTION(id)],SZ_FNAME,"[%d,%d,*]") + } + call pargi (ID_LINE(id,1)) + call pargi (ID_LINE(id,2)) + case 4: + call error (1, "Image dimension greater than 3 not allowed") + } + + call sfree (sp) + } + } + + # Parse the image section. + x1[1] = 1; x2[1] = IM_LEN(im,1); xs[1] = 1 + x1[2] = 1; x2[2] = IM_LEN(im,2); xs[2] = 1 + x1[3] = 1; x2[3] = IM_LEN(im,3); xs[3] = 1 + call id_section (Memc[ID_SECTION(id)], x1, x2, xs, 3) + + # Set the axes. The axis to be identified is the longest one. + i = 1 + if (IM_NDIM(im) > 1 && abs (x1[2]-x2[2]) >= abs (x1[i]-x2[i])) + i = 2 + if (IM_NDIM(im) > 2 && abs (x1[3]-x2[3]) >= abs (x1[i]-x2[i])) + i = 3 + if (IM_NDIM(im) > 3) + call error (1, "Image dimension greater than 3 not allowed") + + switch (i) { + case 1: + j = 2 + k = 3 + case 2: + j = 1 + k = 3 + case 3: + j = 1 + k = 2 + } + + ID_LINE(id,1) = (x1[j] + x2[j]) / 2 + ID_LINE(id,2) = (x1[k] + x2[k]) / 2 + ID_MAXLINE(id,1) = IM_LEN(im, j) + ID_MAXLINE(id,2) = IM_LEN(im, k) + ID_NSUM(id,1) = min (ID_MAXLINE(id,1), ID_NSUM(id,1)) + ID_NSUM(id,2) = min (ID_MAXLINE(id,2), ID_NSUM(id,2)) + call smw_daxis (NULL, NULL, i, ID_NSUM(id,1), ID_NSUM(id,2)) + + call id_maphdr (id, im) + + # Open the image READ_WRITE if possible in order to add REFSPEC. + # This is not done earlier to avoid updating of the WCS. + + call imunmap (im) + if (imaccess (Memc[ID_IMAGE(id)], READ_WRITE) == YES) + im = immap (Memc[ID_IMAGE(id)], READ_WRITE, 0) + else + im = immap (Memc[ID_IMAGE(id)], READ_ONLY, 0) + IM(ID_SH(id)) = im +end + + +# ID_MAPHDR -- Map image header. + +procedure id_maphdr (id, im) + +pointer id # ID pointer +pointer im # IMIO pointer + +int i +pointer mw, sh, smw_openim(), smw_sctran() +errchk smw_openim(), shdr_open(), smw_sctran + +begin + mw = smw_openim (im) + if (SMW_TRANS(mw) == YES) { + if (SMW_PAXIS(mw,1) == 1) + call smw_daxis (mw, im, 2, INDEFI, INDEFI) + else + call smw_daxis (mw, im, 1, INDEFI, INDEFI) + call smw_saxes (mw, NULL, im) + } + call shdr_open (im, mw, ID_LINE(id,1), ID_LINE(id,2), + INDEFI, SHHDR, ID_SH(id)) + sh = ID_SH(id) + + if (SMW_FORMAT(mw) == SMW_MS || SMW_FORMAT(mw) == SMW_ES) { + ID_MAXLINE(id,1) = IM_LEN(im,2) + ID_MAXLINE(id,2) = IM_LEN(im,3) + ID_NSUM(id,1) = 1 + ID_NSUM(id,2) = 1 + ID_LINE(id,1) = max (1, min (ID_MAXLINE(id,1), ID_LINE(id,1))) + ID_LINE(id,2) = 1 + call mfree (ID_APS(id), TY_INT) + call malloc (ID_APS(id), ID_MAXLINE(id,1), TY_INT) + do i = 1, ID_MAXLINE(id,1) { + call shdr_open (im, mw, i, 1, INDEFI, SHHDR, sh) + Memi[ID_APS(id)+i-1] = AP(sh) + } + ID_AP(id,1) = Memi[ID_APS(id)+ID_LINE(id,1)-1] + ID_AP(id,2) = 1 + } else { + call mfree (ID_APS(id), TY_INT) + ID_AP(id,1) = ID_LINE(id,1) + ID_AP(id,2) = ID_LINE(id,2) + } + ID_NPTS(id) = IM_LEN(im, SMW_LAXIS(mw,1)) + + call gt_sets (ID_GT(id), GTXLABEL, LABEL(sh)) + call ic_pstr (ID_IC(id), "ylabel", LABEL(sh)) + call gt_sets (ID_GT(id), GTXUNITS, UNITS(sh)) + call ic_pstr (ID_IC(id), "yunits", UNITS(sh)) + + # Set logical / physical transformations + i = 2 ** (SMW_PAXIS(mw,1) - 1) + ID_LP(id) = smw_sctran (mw, "logical", "physical", i) + ID_PL(id) = smw_sctran (mw, "physical", "logical", i) +end + + +# ID_SECTION -- Parse an image section into its elements. +# 1. The default values must be set by the caller. +# 2. A null image section is OK. +# 3. The first nonwhitespace character must be '['. +# 4. The last interpreted character must be ']'. +# +# This procedure should be replaced with an IMIO procedure at some +# point. + +procedure id_section (section, x1, x2, xs, ndim) + +char section[ARB] # Image section +int x1[ndim] # Starting pixel +int x2[ndim] # Ending pixel +int xs[ndim] # Step +int ndim # Number of dimensions + +int i, ip, a, b, c, temp, ctoi() +define error_ 99 + +begin + # Decode the section string. + ip = 1 + while (IS_WHITE(section[ip])) + ip = ip + 1 + if (section[ip] == '[') + ip = ip + 1 + else if (section[ip] == EOS) + return + else + goto error_ + + do i = 1, ndim { + while (IS_WHITE(section[ip])) + ip = ip + 1 + if (section[ip] == ']') + break + + # Default values + a = x1[i] + b = x2[i] + c = xs[i] + + # Get a:b:c. Allow notation such as "-*:c" + # (or even "-:c") where the step is obviously negative. + + if (ctoi (section, ip, temp) > 0) { # a + a = temp + if (section[ip] == ':') { + ip = ip + 1 + if (ctoi (section, ip, b) == 0) # a:b + goto error_ + } else + b = a + } else if (section[ip] == '-') { # -* + temp = a + a = b + b = temp + ip = ip + 1 + if (section[ip] == '*') + ip = ip + 1 + } else if (section[ip] == '*') # * + ip = ip + 1 + if (section[ip] == ':') { # ..:step + ip = ip + 1 + if (ctoi (section, ip, c) == 0) + goto error_ + else if (c == 0) + goto error_ + } + if (a > b && c > 0) + c = -c + + x1[i] = a + x2[i] = b + xs[i] = c + + while (IS_WHITE(section[ip])) + ip = ip + 1 + if (section[ip] == ',') + ip = ip + 1 + } + + if (section[ip] != ']') + goto error_ + + return +error_ + call error (0, "Error in image section specification") +end diff --git a/noao/rv/rvidlines/idmark.x b/noao/rv/rvidlines/idmark.x new file mode 100644 index 00000000..8f874560 --- /dev/null +++ b/noao/rv/rvidlines/idmark.x @@ -0,0 +1,97 @@ +include <gset.h> +include <smw.h> +include "identify.h" + +procedure id_mark (id, feature) + +pointer id # ID pointer +int feature + +int pix, color, gstati() +real x, y +real mx, my, x1, x2, y1, y2, tick, gap +pointer sp, format, label, ptr +double smw_c1trand() + +define TICK .03 # Tick size in NDC +define GAP .02 # Gap size in NDC + +begin + call ggwind (ID_GP(id), x1, x2, y1, y2) + + x = FIT(id,feature) + + if ((x < min (x1, x2)) || (x > max (x1, x2))) + return + + pix = smw_c1trand (ID_PL(id), PIX(id,feature)) - NP1(ID_SH(id)) + 1 + pix = max (1, min (pix, ID_NPTS(id)-1)) + + call smark (sp) + call salloc (format, SZ_LINE, TY_CHAR) + call salloc (label, SZ_LINE, TY_CHAR) + switch (FTYPE(id,feature)) { + case EMISSION, GEMISSION: + y = max (IMDATA(id,pix), IMDATA(id,pix+1)) + tick = TICK + gap = GAP + call strcpy ("u=180;h=c;v=b;s=0.5", Memc[format], SZ_LINE) + case ABSORPTION, GABSORPTION: + y = min (IMDATA(id,pix), IMDATA(id,pix+1)) + tick = -TICK + gap = -GAP + call strcpy ("u=0;h=c;v=t;s=0.5", Memc[format], SZ_LINE) + } + + call gctran (ID_GP(id), x, y, mx, my, 1, 0) + call gctran (ID_GP(id), mx, my + gap, x1, y1, 0, 1) + call gctran (ID_GP(id), mx, my + gap + tick, x1, y2, 0, 1) + color = gstati (ID_GP(id), G_PLCOLOR) + call gseti (ID_GP(id), G_PLCOLOR, color+1) + call gline (ID_GP(id), x1, y1, x1, y2) + call gseti (ID_GP(id), G_PLCOLOR, color) + + call gctran (ID_GP(id), mx, my + tick + 2 * gap, x1, y2, 0, 1) + color = gstati (ID_GP(id), G_TXCOLOR) + call gseti (ID_GP(id), G_TXCOLOR, color+1) + switch (ID_LABELS(id)) { + case 2: + call sprintf (Memc[label], SZ_LINE, "%d") + call pargi (feature) + call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format]) + case 3: + call sprintf (Memc[label], SZ_LINE, "%0.2f") + call pargd (PIX(id,feature)) + call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format]) + case 4: + if (!IS_INDEFD (USER(id,feature))) { + call sprintf (Memc[label], SZ_LINE, "%0.4f") + call pargd (USER(id,feature)) + call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format]) + } + case 5: + label = Memi[ID_LABEL(id)+feature-1] + if (label != NULL) + call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format]) + case 6: + Memc[label] = EOS + ptr = Memi[ID_LABEL(id)+feature-1] + if (!IS_INDEFD (USER(id,feature))) { + if (ptr != NULL) { + call sprintf (Memc[label], SZ_LINE, "%0.4f %s") + call pargd (USER(id,feature)) + call pargstr (Memc[ptr]) + } else { + call sprintf (Memc[label], SZ_LINE, "%0.4f") + call pargd (USER(id,feature)) + } + } else if (ptr != NULL) + call strcpy (Memc[ptr], Memc[label], SZ_LINE) + if (Memc[label] != EOS) + call gtext (ID_GP(id), x1, y2, Memc[label], Memc[format]) + } + call gseti (ID_GP(id), G_TXCOLOR, color) + + call sfree (sp) + call gflush (ID_GP(id)) +end diff --git a/noao/rv/rvidlines/idnearest.x b/noao/rv/rvidlines/idnearest.x new file mode 100644 index 00000000..41aa4c61 --- /dev/null +++ b/noao/rv/rvidlines/idnearest.x @@ -0,0 +1,29 @@ +include "identify.h" + +# ID_NEAREST -- Find the nearest feature to a given coordinate. + +procedure id_nearest (id, fitnear) + +pointer id # ID pointer +double fitnear # Coordinate to find nearest feature + +int i +double delta, delta1 + +begin + if (ID_NFEATURES(id) < 1) { + ID_CURRENT(id) = 0 + return + } + + ID_CURRENT(id) = 1 + delta = abs (FIT(id,1) - fitnear) + + do i = 2, ID_NFEATURES(id) { + delta1 = abs (FIT(id,i) - fitnear) + if (delta1 < delta) { + ID_CURRENT(id) = i + delta = delta1 + } + } +end diff --git a/noao/rv/rvidlines/idnewfeature.x b/noao/rv/rvidlines/idnewfeature.x new file mode 100644 index 00000000..efa489b4 --- /dev/null +++ b/noao/rv/rvidlines/idnewfeature.x @@ -0,0 +1,87 @@ +include <mach.h> +include "identify.h" + +# ID_NEWFEATURE -- Allocate and initialize memory for a new feature. + +procedure id_newfeature (id, pix, fit, user, wt, width, type, label) + +pointer id # ID pointer +double pix # Pixel coordinate +double fit # Fit coordinate +double user # User coordinate +double wt # Feature weight +real width # Feature width +int type # Feature type +pointer label # Pointer to feature label + +int i, current, strlen() +double delta + +define NALLOC 20 # Length of additional allocations + +begin + if (IS_INDEFD (pix)) + return + + delta = MAX_REAL + do i = 1, ID_NFEATURES(id) { + if (abs (pix - PIX(id,i)) < delta) { + delta = abs (pix - PIX(id,i)) + current = i + } + } + + if (delta >= ID_MINSEP(id)) { + ID_NFEATURES(id) = ID_NFEATURES(id) + 1 + if (ID_NALLOC(id) < ID_NFEATURES(id)) { + ID_NALLOC(id) = ID_NALLOC(id) + NALLOC + call realloc (ID_PIX(id), ID_NALLOC(id), TY_DOUBLE) + call realloc (ID_FIT(id), ID_NALLOC(id), TY_DOUBLE) + call realloc (ID_USER(id), ID_NALLOC(id), TY_DOUBLE) + call realloc (ID_WTS(id), ID_NALLOC(id), TY_DOUBLE) + call realloc (ID_FWIDTHS(id), ID_NALLOC(id), TY_REAL) + call realloc (ID_FTYPES(id), ID_NALLOC(id), TY_INT) + call realloc (ID_LABEL(id), ID_NALLOC(id), TY_POINTER) + call aclri (Memi[ID_LABEL(id)+ID_NALLOC(id)-NALLOC], NALLOC) + } + for (current=ID_NFEATURES(id); (current>1)&&(pix<PIX(id,current-1)); + current=current-1) { + PIX(id,current) = PIX(id,current-1) + FIT(id,current) = FIT(id,current-1) + USER(id,current) = USER(id,current-1) + WTS(id,current) = WTS(id,current-1) + FWIDTH(id,current) = FWIDTH(id,current-1) + FTYPE(id,current) = FTYPE(id,current-1) + Memi[ID_LABEL(id)+current-1] = Memi[ID_LABEL(id)+current-2] + } + PIX(id,current) = pix + FIT(id,current) = fit + USER(id,current) = user + WTS(id,current) = wt + FWIDTH(id,current) = width + FTYPE(id,current) = type + if (label != NULL) { + i = strlen (Memc[label]) + call malloc (Memi[ID_LABEL(id)+current-1], i, TY_CHAR) + call strcpy (Memc[label], Memc[Memi[ID_LABEL(id)+current-1]], i) + } else + Memi[ID_LABEL(id)+current-1] = NULL + ID_NEWFEATURES(id) = YES + } else if (abs (fit-user) < abs (FIT(id,current)-USER(id,current))) { + PIX(id,current) = pix + FIT(id,current) = fit + USER(id,current) = user + WTS(id,current) = wt + FWIDTH(id,current) = width + FTYPE(id,current) = type + if (label != NULL) { + i = strlen (Memc[label]) + call malloc (Memi[ID_LABEL(id)+current-1], i, TY_CHAR) + call strcpy (Memc[label], Memc[Memi[ID_LABEL(id)+current-1]], i) + } else + Memi[ID_LABEL(id)+current-1] = NULL + ID_NEWFEATURES(id) = YES + } + + ID_CURRENT(id) = current +end diff --git a/noao/rv/rvidlines/idnoextn.x b/noao/rv/rvidlines/idnoextn.x new file mode 100644 index 00000000..6c82d778 --- /dev/null +++ b/noao/rv/rvidlines/idnoextn.x @@ -0,0 +1,11 @@ +# ID_NOEXTN -- Remove standard image extensions. + +procedure id_noextn (image) + +char image[ARB] # Image name + +int strlen() + +begin + call xt_imroot (image, image, strlen (image)) +end diff --git a/noao/rv/rvidlines/idpeak.x b/noao/rv/rvidlines/idpeak.x new file mode 100644 index 00000000..9cba49c4 --- /dev/null +++ b/noao/rv/rvidlines/idpeak.x @@ -0,0 +1,23 @@ +include "identify.h" + +# ID_PEAK -- Find the peak value above continuum. + +double procedure id_peak (id, pix) + +pointer id # ID pointer +double pix # Pixel position +double peak # Peak value + +int c, l, u + +begin + if (IS_INDEFD(pix)) + return (INDEFD) + + c = nint (pix) + l = max (1, nint (pix - ID_FWIDTH(id))) + u = min (ID_NPTS(id), nint (pix + ID_FWIDTH(id))) + peak = IMDATA(id,c) - (IMDATA(id,l) + IMDATA(id,u)) / 2. + + return (peak) +end diff --git a/noao/rv/rvidlines/idrms.x b/noao/rv/rvidlines/idrms.x new file mode 100644 index 00000000..257675de --- /dev/null +++ b/noao/rv/rvidlines/idrms.x @@ -0,0 +1,28 @@ +include "identify.h" + +# ID_RMS -- Compute RMS of fit about the user coordinates + +double procedure id_rms (id) + +pointer id # ID pointer + +int i, nrms +double rms, id_zshiftd() + +begin + rms = 0. + nrms = 0 + for (i=1; i<=ID_NFEATURES(id); i=i+1) { + if (!IS_INDEFD (USER(id,i)) && WTS(id,i) != 0.) { + rms = rms + (id_zshiftd (id, FIT(id,i), 0) - USER(id,i)) ** 2 + nrms = nrms + 1 + } + } + + if (nrms > 0) + rms = sqrt (rms / nrms) + else + rms = INDEFD + + return (rms) +end diff --git a/noao/rv/rvidlines/idshift.x b/noao/rv/rvidlines/idshift.x new file mode 100644 index 00000000..3fa9f9bf --- /dev/null +++ b/noao/rv/rvidlines/idshift.x @@ -0,0 +1,65 @@ +include "identify.h" + +define NBIN 10 # Bin parameter for mode determination + +# ID_SHIFT -- Determine a shift by correlating feature user positions +# with peaks in the image data. + +double procedure id_shift (id) + +pointer id # ID pointer + +int i, j, npeaks, ndiff, find_peaks() +real d, dmin +double pix, id_center(), id_fitpt() +pointer x, y, diff +errchk malloc, find_peaks + +begin + # Find the peaks in the image data and center. + call malloc (x, ID_NPTS(id), TY_REAL) + npeaks = find_peaks (IMDATA(id,1), Memr[x], ID_NPTS(id), 0., + int (ID_MINSEP(id)), 0, ID_MAXFEATURES(id), 0., false) + + # Center the peaks and convert to user coordinates. + call malloc (y, npeaks, TY_DOUBLE) + j = 0 + do i = 1, npeaks { + pix = id_center (id, double(Memr[x+i-1]), 1, ID_FWIDTH(id), + ID_FTYPE(id), NO) + if (!IS_INDEFD (pix)) { + Memd[y+j] = id_fitpt (id, pix) + j = j + 1 + } + } + npeaks = j + + # Compute differences with feature list. + ndiff = npeaks * ID_NFEATURES(id) + call malloc (diff, ndiff, TY_REAL) + ndiff = 0 + do i = 1, ID_NFEATURES(id) { + do j = 1, npeaks { + Memr[diff+ndiff] = Memd[y+j-1] - FIT(id,i) + ndiff = ndiff + 1 + } + } + call mfree (x, TY_REAL) + call mfree (y, TY_DOUBLE) + + # Sort the differences and find the mode. + call asrtr (Memr[diff], Memr[diff], ndiff) + + dmin = Memr[diff+ndiff-1] - Memr[diff] + do i = 0, ndiff-NBIN-1 { + j = i + NBIN + d = Memr[diff+j] - Memr[diff+i] + if (d < dmin) { + dmin = d + pix = Memr[diff+i] + d / 2. + } + } + call mfree (diff, TY_REAL) + + return (pix) +end diff --git a/noao/rv/rvidlines/idshow.x b/noao/rv/rvidlines/idshow.x new file mode 100644 index 00000000..a01fba6f --- /dev/null +++ b/noao/rv/rvidlines/idshow.x @@ -0,0 +1,83 @@ +include "identify.h" + +# ID_SHOW -- Show parameter information. + +procedure id_show (id, file) + +pointer id # ID pointer +char file[ARB] # File + +char line[SZ_LINE] +int fd + +int open(), ic_geti() +errchk open() + +begin + fd = open (file, APPEND, TEXT_FILE) + + call sysid (line, SZ_LINE) + call fprintf (fd, "%s\n") + call pargstr (line) + + call fprintf (fd, "image %s\n") + call pargstr (Memc[ID_IMAGE(id)]) + call fprintf (fd, "nsum %d\n") + call pargi (ID_NSUM(id,1)) + switch (ID_FTYPE(id)) { + case EMISSION: + call fprintf (fd, "ftype emission\n") + case ABSORPTION: + call fprintf (fd, "ftype absorption\n") + case GEMISSION: + call fprintf (fd, "ftype gemission\n") + case GABSORPTION: + call fprintf (fd, "ftype gabsorption\n") + } + switch (ID_LABELS(id)) { + case 2: + call fprintf (fd, "labels index\n") + case 3: + call fprintf (fd, "labels pixel\n") + case 4: + call fprintf (fd, "labels coords\n") + case 5: + call fprintf (fd, "labels user\n") + case 6: + call fprintf (fd, "labels both\n") + default: + call fprintf (fd, "labels none\n") + } + call fprintf (fd, "maxfeatures %d\n") + call pargi (ID_MAXFEATURES(id)) + call fprintf (fd, "match %g\n") + call pargr (ID_MATCH(id)) + call fprintf (fd, "zwidth %g\n") + call pargr (ID_ZWIDTH(id)) + call fprintf (fd, "fwidth %g\n") + call pargr (ID_FWIDTH(id)) + call fprintf (fd, "database %s\n") + call pargstr (Memc[ID_DATABASE(id)]) + call fprintf (fd, "coordlist %s\n") + call pargstr (Memc[ID_COORDLIST(id)]) + call fprintf (fd, "cradius %g\n") + call pargr (ID_CRADIUS(id)) + call fprintf (fd, "threshold %g\n") + call pargr (ID_THRESHOLD(id)) + call fprintf (fd, "minsep %g\n") + call pargr (ID_MINSEP(id)) + if (ID_CV(id) != NULL) { + call fprintf (fd, "function = %s\n") + call ic_gstr (ID_IC(id), "function", line, SZ_LINE) + call pargstr (line) + call fprintf (fd, "order = %d\n") + call pargi (ic_geti (ID_IC(id), "order")) + call fprintf (fd, "Fit at first pixel = %0.8g\n") + call pargd (FITDATA(id,1)) + call fprintf (fd, "Average fit interval = %0.8g\n") + call pargd ((FITDATA(id,ID_NPTS(id))-FITDATA(id,1))/ + (ID_NPTS(id)-1)) + } + + call close (fd) +end diff --git a/noao/rv/rvidlines/idvelocity.x b/noao/rv/rvidlines/idvelocity.x new file mode 100644 index 00000000..c62e5c89 --- /dev/null +++ b/noao/rv/rvidlines/idvelocity.x @@ -0,0 +1,188 @@ +include <smw.h> +include <units.h> +include "identify.h" + + +# ID_VELOCITY -- Compute velocity. + +procedure id_velocity (id, interactive) + +pointer id # ID pointer +int interactive # Called interactively? + +int i, n +double z, sumz, sumz2, sumw, zerr, zhelio, v, verr, id_zval() + +begin + sumz = 0 + sumw = 0 + n = 0 + for (i=1; i <= ID_NFEATURES(id); i = i + 1) { + if (IS_INDEFD (USER(id,i)) || WTS(id,i) == 0.) + next + z = id_zval (id, FIT(id,i), USER(id,i)) + sumz = sumz + WTS(id,i) * z + sumw = sumw + WTS(id,i) + n = n + 1 + } + + if (sumw > 0.) { + zhelio = ID_ZHELIO(id) + sumz = sumz / sumw + + sumz2 = 0. + for (i=1; i <= ID_NFEATURES(id); i = i + 1) { + if (IS_INDEFD (USER(id,i)) || WTS(id,i) == 0.) + next + z = id_zval (id, FIT(id,i), USER(id,i)) + sumz2 = sumz2 + WTS(id,i) * (z - sumz) ** 2 + } + if (sumz2 > 0.) + sumz2 = sqrt (sumz2 / sumw) + else + sumz2 = 0. + zerr = sumz2 + if (n > 1) + zerr = zerr / sqrt (n - 1.) + + if (interactive == YES) { + v = (sumz + zhelio) * VLIGHT + verr = zerr * VLIGHT + + if (zhelio == 0D0) + call printf ( + "%s%s: Lines=%3d, Vobs=%.5g (%.5g), Zobs=%.5g (%.5g)\n") + else + call printf ( + "%s%s: Lines=%3d, Vhelio=%.5g (%.5g), Zhelio=%.5g (%.5g)\n") + + call pargstr (Memc[ID_IMAGE(id)]) + call pargstr (Memc[ID_SECTION(id)]) + call pargi (n) + call pargd (v) + call pargd (verr) + call pargd (sumz + zhelio) + call pargd (zerr) + } + ID_REDSHIFT(id) = sumz + ID_RMSRED(id) = sumz2 + ID_ZHELIO(id) = zhelio + } +end + + +# ID_ZVAL -- Compute Z value. + +double procedure id_zval (id, x, xref) + +pointer id #I Identify pointer +double x #I Coordinate +double xref #I Reference coordinate +double z #O Z value + +double y, yref +pointer un + +begin + y = x + yref = xref + + un = UN(ID_SH(id)) + if (UN_LOG(un) == YES) { + y = 10D0 ** y + yref = 10D0 ** yref + } + if (UN_INV(un) == YES) { + y = 1D0 / y + yref = 1D0 / yref + } + + switch (UN_CLASS(un)) { + case UN_WAVE: + z = (y - yref) / yref + case UN_FREQ, UN_ENERGY: + z = (yref - y) / y + case UN_VEL: + y = sqrt ((1 + y) / (1 - y)) + yref = sqrt ((1 + yref) / (1 - yref)) + z = (y - yref) / yref + case UN_DOP: + y = y + 1 + yref = yref + 1 + z = (y - yref) / yref + } + + return (z) +end + + + +# ID_ZSHIFTD -- Shift coordinate by redshift. + +double procedure id_zshiftd (id, x, dir) + +pointer id #I Identify pointer +double x #I Coordinate +int dir #I Direction (0=to rest, 1=from rest) +double y #O Shifted coordinate + +pointer un + +begin + y = x + + un = UN(ID_SH(id)) + if (UN_LOG(un) == YES) + y = 10D0 ** y + if (UN_INV(un) == YES) + y = 1D0 / y + + switch (UN_CLASS(un)) { + case UN_WAVE: + if (dir == 0) + y = y / (1 + ID_REDSHIFT(id)) + else + y = y * (1 + ID_REDSHIFT(id)) + case UN_FREQ, UN_ENERGY: + if (dir == 0) + y = y * (1 + ID_REDSHIFT(id)) + else + y = y / (1 + ID_REDSHIFT(id)) + case UN_VEL: + y = sqrt ((1 + y) / (1 - y)) + if (dir == 0) + y = y / (1 + ID_REDSHIFT(id)) + else + y = y * (1 + ID_REDSHIFT(id)) + y = y ** 2 + y = (y - 1) / (y + 1) + case UN_DOP: + y = (y + 1) + if (dir == 0) + y = (y + 1) / (1 + ID_REDSHIFT(id)) - 1 + else + y = (y + 1) * (1 + ID_REDSHIFT(id)) - 1 + } + + if (UN_INV(un) == YES) + y = 1D0 / y + if (UN_LOG(un) == YES) + y = log10 (y) + + return (y) +end + + +# ID_ZSHIFTR -- Shift coordinate by redshift. + +real procedure id_zshiftr (id, x, dir) + +pointer id #I Identify pointer +real x #I Coordinate +int dir #I Direction (0=to rest, 1=from rest) + +double id_zshiftd() + +begin + return (real (id_zshiftd (id, double(x), dir))) +end diff --git a/noao/rv/rvidlines/idvhelio.x b/noao/rv/rvidlines/idvhelio.x new file mode 100644 index 00000000..803efc2d --- /dev/null +++ b/noao/rv/rvidlines/idvhelio.x @@ -0,0 +1,102 @@ +include <error.h> + +# ID_VHELIO -- Compute helocentric velocity. + +procedure id_vhelio (im, vhelio, hjd, fd) + +pointer im #I IMIO pointer +double vhelio #O Heliocentric velocity correction +double hjd #O Heliocentric Julian Date +int fd #I Log file descriptor + +bool newobs, obshead +int year, month, day, flags, dummy +double ra, dec, ep, ut, lt +double epoch, vrot, vbary, vorb +double latitude, longitude, altitude +pointer sp, str1, str2, tmp, obs, kp + +int dtm_decode() +double imgetd(), obsgetd() +pointer clopset() + +int err +data err/0/ + +errchk imgetd, imgstr, obsopen + + +begin + call smark (sp) + call salloc (str1, SZ_LINE, TY_CHAR) + call salloc (str2, SZ_LINE, TY_CHAR) + + obs = NULL + + iferr { + # Get the observatory data. + call clgstr ("observatory", Memc[str1], SZ_LINE) + tmp = NULL + call obsimopen (tmp, im, Memc[str1], NO, newobs, obshead) + obs = tmp + + latitude = obsgetd (obs, "latitude") + longitude = obsgetd (obs, "longitude") + altitude = obsgetd (obs, "altitude") + + # Get the image header data. + kp = clopset ("keywpars") + + call clgpset (kp, "date_obs", Memc[str1], SZ_LINE) + call imgstr (im, Memc[str1], Memc[str2], SZ_LINE) + if (dtm_decode (Memc[str2],year,month,day,ut,flags) == ERR) + call error (1, "Error in date string") + + call clgpset (kp, "ut", Memc[str1], SZ_LINE) + call imgstr (im, Memc[str1], Memc[str2], SZ_LINE) + if (dtm_decode (Memc[str2],dummy,dummy,dummy,ut,flags) == ERR) { + iferr (ut = imgetd (im, Memc[str1])) + call error (1, "Error in UT keyword") + } + + call clgpset (kp, "ra", Memc[str1], SZ_LINE) + ra = imgetd (im, Memc[str1]) + call clgpset (kp, "dec", Memc[str1], SZ_LINE) + dec = imgetd (im, Memc[str1]) + call clgpset (kp, "epoch", Memc[str1], SZ_LINE) + ep = imgetd (im, Memc[str1]) + call clcpset (kp) + + # Determine epoch of observation and precess coordinates. + call ast_date_to_epoch (year, month, day, ut, epoch) + call ast_precess (ra, dec, ep, ra, dec, epoch) + + # Determine velocity components. + call ast_vorbit (ra, dec, epoch, vorb) + call ast_vbary (ra, dec, epoch, vbary) + call ast_vrotate (ra, dec, epoch, latitude, longitude, + altitude, vrot) + call ast_hjd (ra, dec, epoch, lt, hjd) + + vhelio = vrot + vbary + vorb + + if (fd != NULL) + call obslog (obs, + "RVIDLINES", "latitude longitude altitude", fd) + + } then { + vhelio = 0. + if (err == 0) { + call eprintf ("Warning: Can't compute heliocentric velocity\n") + call erract (EA_WARN) + err = err + 1 + } + } + + iferr { # This IFERR is to clear errcode. + if (obs != NULL) + call obsclose (obs) + } then + ; + call sfree (sp) +end diff --git a/noao/rv/rvidlines/mkpkg b/noao/rv/rvidlines/mkpkg new file mode 100644 index 00000000..a18105f8 --- /dev/null +++ b/noao/rv/rvidlines/mkpkg @@ -0,0 +1,47 @@ +# RVIDLINES + +$checkout libpkg.a .. +$update libpkg.a +$checkin libpkg.a .. +$exit + +libpkg.a: + $ifeq (USE_GENERIC, yes) + $ifolder (peaks.x, peaks.gx) + $generic -k peaks.gx -o peaks.x $endif $endif + + # Version from V2.10.3 ONEDSPEC that conflicts with V2.10.2 RV + iddeblend.x <mach.h> + + idcenter.x identify.h <gset.h> <smw.h> + idcolon.x identify.h <error.h> <gset.h> <smw.h> + iddb.x identify.h <imset.h> <math/curfit.h> <smw.h> + iddelete.x identify.h + iddofit.x identify.h + iddoshift.x identify.h + idfitdata.x identify.h <smw.h> + idfixx.x <smw.h> + idgdata.x identify.h <imhdr.h> <imio.h> <pkg/gtools.h> <smw.h> + idgraph.x identify.h <gset.h> <pkg/gtools.h> <smw.h> + ididentify.x identify.h <error.h> <gset.h> <imhdr.h> <smw.h> + idinit.x identify.h <gset.h> <math/curfit.h> + idlabel.x + idlinelist.x identify.h <error.h> <mach.h> + idlog.x identify.h <smw.h> <time.h> + idmap.x identify.h <ctype.h> <imhdr.h> <pkg/gtools.h> <smw.h> + idmark.x identify.h <gset.h> <smw.h> + idnearest.x identify.h + idnewfeature.x identify.h <mach.h> + idnoextn.x + idpeak.x identify.h + idrms.x identify.h + idshift.x identify.h + idshow.x identify.h + idvelocity.x identify.h <smw.h> <units.h> + idvhelio.x <error.h> + peaks.x + reidentify.x identify.h <error.h> <gset.h> <imhdr.h> + t_identify.x identify.h <mach.h> <pkg/gtools.h> + t_reidentify.x identify.h <error.h> <fset.h> <gset.h> <pkg/gtools.h>\ + <smw.h> + ; diff --git a/noao/rv/rvidlines/peaks.gx b/noao/rv/rvidlines/peaks.gx new file mode 100644 index 00000000..32764505 --- /dev/null +++ b/noao/rv/rvidlines/peaks.gx @@ -0,0 +1,447 @@ +# PEAKS -- The following procedures are general numerical functions +# dealing with finding peaks in a data array. +# +# FIND_PEAKS Find the peaks in the data array. +# FIND_LOCAL_MAXIMA Find the local minima/maxima in the data array. +# IS_LOCAL_MAX Test a point to determine if it is a local min/max. +# FIND_THRESHOLD Find the peaks with positions satisfying threshold +# and contrast constraints. +# FIND_ISOLATED Flag peaks which are within separation of a peak +# with a higher peak value. +# FIND_NMAX Select up to the nmax highest ranked peaks. +# COMPARE Compare procedure for sort used in FIND_PEAKS. + +# FIND_PEAKS -- Find the peaks in the data array. +# +# The peaks are found using the following algorithm: +# +# 1. Find the local minima/maxima. +# 2. Reject peaks below the threshold. +# 3. Determine the ranks of the remaining peaks. +# 4. Flag weaker peaks within separation of a stronger peak. +# 5. Accept at most the nmax strongest peaks. +# +# Indefinite points are ignored. The peak positions are returned in the +# array x. + +$for (r) +int procedure find_peaks (data, x, npoints, contrast, separation, edge, nmax, + threshold, debug) + +# Procedure parameters: +PIXEL data[npoints] # Input data array +PIXEL x[npoints] # Output peak position array +int npoints # Number of data points +real contrast # Maximum contrast between strongest and weakest +int separation # Minimum separation between peaks +int edge # Minimum distance from the edge +int nmax # Maximum number of peaks to be returned +real threshold # Minimum threshold level for peaks +bool debug # Print diagnostic information? + +int i, j +int nlmax, nthreshold, nisolated, npeaks +pointer sp, rank + +int find_local_maxima(), find_threshold(), find_isolated(), find_nmax() +int compare() + +extern compare() + +pointer y +int maxima +common /sort/ y, maxima + +begin + # Find the local minima/maxima in data and put column positions in x.. + if (contrast < 0.) + maxima = NO + else + maxima = YES + + nlmax = find_local_maxima (data, x, npoints, debug) + + # Reject local minima/maxima near the edge. + if (edge > 0) { + j = 0 + do i = 1, nlmax { + if ((x[i] > edge) && (x[i] <= npoints - edge)) { + j = j + 1 + x[j] = x[i] + } + } + nlmax = j + } + + # Allocate a working array y. + call smark (sp) + call salloc (y, npoints, TY_PIXEL) + + # Reject the local minima/maxima which do not satisfy the thresholds. + # The array y is set to the peak values of the remaining peaks. + nthreshold = find_threshold (data, x, Mem$t[y], nlmax, + contrast, threshold, debug) + + # Rank the peaks by peak value. + call salloc (rank, nthreshold, TY_INT) + do i = 1, nthreshold + Memi[rank + i - 1] = i + call qsort (Memi[rank], nthreshold, compare) + + # Reject the weaker peaks within sep of a stronger peak. + nisolated = find_isolated (x, Memi[rank], nthreshold, separation, + debug) + + # Select the strongest nmax peaks. + npeaks = find_nmax (data, x, Memi[rank], nthreshold, nmax, debug) + + call sfree (sp) + return (npeaks) +end + + +# FIND_LOCAL_MAXIMA -- Find the local minima/maxima in the data array. +# +# A data array is input and the local minima/maxima positions array is output. +# The number of local minima/maxima found is returned. + +int procedure find_local_maxima (data, x, npoints, debug) + +PIXEL data[npoints] # Input data array +PIXEL x[npoints] # Output local min/max positions array +int npoints # Number of input points +bool debug # Print debugging information? + +int i, nlmax + +bool is_local_max() + +begin + nlmax = 0 + do i = 1, npoints { + if (is_local_max (i, data, npoints)) { + nlmax = nlmax + 1 + x[nlmax] = i + } + } + + if (debug) { + call printf (" Number of local minima/maxima found = %d.\n") + call pargi (nlmax) + } + + return (nlmax) +end + + +# IS_LOCAL_MAX -- Test a point to determine if it is a local minima/maximum. +# +# Indefinite points are ignored. + +bool procedure is_local_max (index, data, npoints) + +# Procedure parameters: +int index # Index to test for local minimum/maximum +PIXEL data[npoints] # Data values +int npoints # Number of points in the data vector + +int i, j, nright, nleft + +pointer y +int maxima +common /sort/ y, maxima + +begin + # INDEF points cannot be local minima/axima. + if (IS_INDEF (data[index])) + return (FALSE) + + # Find the left and right indices where data values change and the + # number of points with the same value. Ignore INDEF points. + nleft = 0 + for (i = index - 1; i >= 1; i = i - 1) { + if (!IS_INDEF (data[i])) { + if (data[i] != data[index]) + break + nleft = nleft + 1 + } + } + nright = 0 + for (j = index + 1; i <= npoints; j = j + 1) { + if (!IS_INDEF (data[j])) { + if (data[j] != data[index]) + break + nright = nright + 1 + } + } + + # Test for failure to be a local minima/axima + if (maxima == YES) { + if ((i == 0) && (j == npoints)) { + return (FALSE) # Data is constant + } else if (i == 0) { + if (data[j] > data[index]) + return (FALSE) # Data increases to right + } else if (j == npoints) { + if (data[i] > data[index]) # Data increase to left + return (FALSE) + } else if ((data[i] > data[index]) || (data[j] > data[index])) { + return (FALSE) # Not a local maximum + } else if (!((nleft - nright == 0) || (nleft - nright == 1))) { + return (FALSE) # Not center of plateau + } + } else { + if ((i == 0) && (j == npoints)) { + return (FALSE) # Data is constant + } else if (i == 0) { + if (data[j] < data[index]) + return (FALSE) # Data decreases to right + } else if (j == npoints) { + if (data[i] < data[index]) # Data decrease to left + return (FALSE) + } else if ((data[i] < data[index]) || (data[j] < data[index])) { + return (FALSE) # Not a local maximum + } else if (!((nleft - nright == 0) || (nleft - nright == 1))) { + return (FALSE) # Not center of plateau + } + } + + # Point is a local minima/maxima + return (TRUE) +end + + +# FIND_THRESHOLD -- Find the peaks with positions satisfying threshold +# and contrast constraints. +# +# The input is the data array, data, and the peak positions array, x. +# The x array is resorted to the nthreshold peaks satisfying the constraints. +# The corresponding nthreshold data values are returned the y array. +# The number of peaks satisfying the constraints (nthreshold) is returned. + +int procedure find_threshold (data, x, y, npoints, contrast, threshold, debug) + +PIXEL data[ARB] # Input data values +PIXEL x[npoints] # Input/Output peak positions +PIXEL y[npoints] # Output peak data values +int npoints # Number of peaks input +real contrast # Contrast constraint +real threshold # Threshold constraint +bool debug # Print debugging information? + +int i, j, nthreshold +PIXEL minval, maxval, lcut + +pointer dummy +int maxima +common /sort/ dummy, maxima + +begin + # Set the y array to be the values at the peak positions. + do i = 1, npoints { + j = x[i] + y[i] = data[j] + } + + # Determine the min and max values of the peaks. + call alim$t (y, npoints, minval, maxval) + + # Set the threshold based on the max of the absolute threshold and the + # contrast. Use arlt to set peaks below threshold to INDEF. + + if (maxima == YES) { + lcut = max (PIXEL (threshold), PIXEL (contrast * maxval)) + call arlt$t (y, npoints, lcut, INDEF) + } else { + lcut = threshold + call argt$t (y, npoints, lcut, INDEF) + } + + if (debug) { + call printf (" Highest peak value = %g.\n") + call parg$t (maxval) + call printf (" Peak cutoff threshold = %g.\n") + call parg$t (lcut) + do i = 1, npoints { + if (IS_INDEF (y[i])) { + j = x[i] + call printf ( + " Peak at column %d with value %g below threshold.\n") + call pargi (j) + call parg$t (data[j]) + } + } + } + + # Determine the number of acceptable peaks & resort the x and y arrays. + nthreshold = 0 + do i = 1, npoints { + if (IS_INDEF (y[i])) + next + nthreshold = nthreshold + 1 + x[nthreshold] = x[i] + y[nthreshold] = y[i] + } + + if (debug) { + call printf (" Number of peaks above the threshold = %d.\n") + call pargi (nthreshold) + } + + return (nthreshold) +end + +# FIND_ISOLATED -- Flag peaks which are within separation of a peak +# with a higher peak value. +# +# The peak positions, x, and their ranks, rank, are input. +# The rank array contains the indices of the peak positions in order from +# the highest peak value to the lowest peak value. Starting with +# highest rank (rank[1]) all peaks of lower rank within separation +# are marked by setting their positions to INDEF. The number of +# unflaged peaks is returned. + +int procedure find_isolated (x, rank, npoints, separation, debug) + +# Procedure parameters: +PIXEL x[npoints] # Positions of points +int rank[npoints] # Rank of peaks +int npoints # Number of peaks +int separation # Minimum allowed separation +bool debug # Print diagnostic information + +int i, j +int nisolated + +begin + # Eliminate close neighbors. The eliminated + # peaks are marked by setting their positions to INDEF. + nisolated = 0 + do i = 1, npoints { + if (IS_INDEF (x[rank[i]])) + next + nisolated = nisolated + 1 + do j = i + 1, npoints { + if (IS_INDEF (x[rank[j]])) + next + if (abs (x[rank[i]] - x[rank[j]]) < separation) { + if (debug) { + call printf ( + " Peak at column %d too near peak at column %d.\n") + call pargi (int (x[rank[j]])) + call pargi (int (x[rank[i]])) + } + x[rank[j]] = INDEF + } + } + } + + if (debug) { + call printf (" Number of peaks separated by %d pixels = %d.\n") + call pargi (separation) + call pargi (nisolated) + } + + # Return number of isolated peaks. + return (nisolated) +end + + +# FIND_NMAX -- Select up to the nmax highest ranked peaks. +# +# The data values, data, peak positions, x, and their ranks, rank, are input. +# The data values are used only in printing debugging information. +# Peak positions previously eliminated are flaged by the value INDEF. +# The rank array contains the indices to the peak positions in order from +# the highest peak value to the lowest peak value. +# First all but the nmax highest ranked peaks (which have not been previously +# eliminated) are eliminated by marking their positions with the value INDEF. +# Then the remaining peaks are resorted to contain only the unflaged +# peaks and the number of such peaks is returned. + +int procedure find_nmax (data, x, rank, npoints, nmax, debug) + +PIXEL data[ARB] # Input data values +PIXEL x[npoints] # Peak positions +int rank[npoints] # Ranks of peaks +int npoints # Number of input peaks +int nmax # Max number of peaks to be selected +bool debug # Print debugging information? + +int i, j, npeaks + +begin + # Only mark peaks to reject if the number peaks is greater than nmax. + if (nmax < npoints) { + npeaks = 0 + do i = 1, npoints { + if (IS_INDEF (x[rank[i]])) + next + npeaks = npeaks + 1 + if (npeaks > nmax) { + if (debug) { + j = x[rank[i]] + call printf ( + " Reject peak at column %d with rank %d and value %g.\n") + call pargi (j) + call pargi (i) + call parg$t (data[j]) + } + x[rank[i]] = INDEF + } + } + } + + # Eliminate INDEF points and determine the number of spectra found. + npeaks = 0 + do i = 1, npoints { + if (IS_INDEF (x[i])) + next + npeaks = npeaks + 1 + x[npeaks] = x[i] + } + + return (npeaks) +end + + +# COMPARE -- Compare procedure for sort used in FIND_PEAKS. +# Larger values are indexed first. INDEF values are indexed last. + +int procedure compare (index1, index2) + +# Procedure parameters: +int index1 # Comparison index +int index2 # Comparison index + +pointer y +int maxima +common /sort/ y, maxima + +begin + # INDEF points are considered to be smallest/largest possible values. + if (maxima == YES) { + if (IS_INDEF (Mem$t[y - 1 + index1])) + return (1) + else if (IS_INDEF (Mem$t[y - 1 + index2])) + return (-1) + else if (Mem$t[y - 1 + index1] < Mem$t[y - 1 + index2]) + return (1) + else if (Mem$t[y - 1 + index1] > Mem$t[y - 1 + index2]) + return (-1) + else + return (0) + } else { + if (IS_INDEF (Mem$t[y - 1 + index1])) + return (-1) + else if (IS_INDEF (Mem$t[y - 1 + index2])) + return (1) + else if (Mem$t[y - 1 + index1] < Mem$t[y - 1 + index2]) + return (-1) + else if (Mem$t[y - 1 + index1] > Mem$t[y - 1 + index2]) + return (1) + else + return (0) + } +end +$endfor diff --git a/noao/rv/rvidlines/peaks.x b/noao/rv/rvidlines/peaks.x new file mode 100644 index 00000000..8ee27f51 --- /dev/null +++ b/noao/rv/rvidlines/peaks.x @@ -0,0 +1,446 @@ +# PEAKS -- The following procedures are general numerical functions +# dealing with finding peaks in a data array. +# +# FIND_PEAKS Find the peaks in the data array. +# FIND_LOCAL_MAXIMA Find the local minima/maxima in the data array. +# IS_LOCAL_MAX Test a point to determine if it is a local min/max. +# FIND_THRESHOLD Find the peaks with positions satisfying threshold +# and contrast constraints. +# FIND_ISOLATED Flag peaks which are within separation of a peak +# with a higher peak value. +# FIND_NMAX Select up to the nmax highest ranked peaks. +# COMPARE Compare procedure for sort used in FIND_PEAKS. + +# FIND_PEAKS -- Find the peaks in the data array. +# +# The peaks are found using the following algorithm: +# +# 1. Find the local minima/maxima. +# 2. Reject peaks below the threshold. +# 3. Determine the ranks of the remaining peaks. +# 4. Flag weaker peaks within separation of a stronger peak. +# 5. Accept at most the nmax strongest peaks. +# +# Indefinite points are ignored. The peak positions are returned in the +# array x. + + +int procedure find_peaks (data, x, npoints, contrast, separation, edge, nmax, + threshold, debug) + +# Procedure parameters: +real data[npoints] # Input data array +real x[npoints] # Output peak position array +int npoints # Number of data points +real contrast # Maximum contrast between strongest and weakest +int separation # Minimum separation between peaks +int edge # Minimum distance from the edge +int nmax # Maximum number of peaks to be returned +real threshold # Minimum threshold level for peaks +bool debug # Print diagnostic information? + +int i, j +int nlmax, nthreshold, nisolated, npeaks +pointer sp, rank + +int find_local_maxima(), find_threshold(), find_isolated(), find_nmax() +int compare() + +extern compare() + +pointer y +int maxima +common /sort/ y, maxima + +begin + # Find the local minima/maxima in data and put column positions in x.. + if (contrast < 0.) + maxima = NO + else + maxima = YES + + nlmax = find_local_maxima (data, x, npoints, debug) + + # Reject local minima/maxima near the edge. + if (edge > 0) { + j = 0 + do i = 1, nlmax { + if ((x[i] > edge) && (x[i] <= npoints - edge)) { + j = j + 1 + x[j] = x[i] + } + } + nlmax = j + } + + # Allocate a working array y. + call smark (sp) + call salloc (y, npoints, TY_REAL) + + # Reject the local minima/maxima which do not satisfy the thresholds. + # The array y is set to the peak values of the remaining peaks. + nthreshold = find_threshold (data, x, Memr[y], nlmax, + contrast, threshold, debug) + + # Rank the peaks by peak value. + call salloc (rank, nthreshold, TY_INT) + do i = 1, nthreshold + Memi[rank + i - 1] = i + call qsort (Memi[rank], nthreshold, compare) + + # Reject the weaker peaks within sep of a stronger peak. + nisolated = find_isolated (x, Memi[rank], nthreshold, separation, + debug) + + # Select the strongest nmax peaks. + npeaks = find_nmax (data, x, Memi[rank], nthreshold, nmax, debug) + + call sfree (sp) + return (npeaks) +end + + +# FIND_LOCAL_MAXIMA -- Find the local minima/maxima in the data array. +# +# A data array is input and the local minima/maxima positions array is output. +# The number of local minima/maxima found is returned. + +int procedure find_local_maxima (data, x, npoints, debug) + +real data[npoints] # Input data array +real x[npoints] # Output local min/max positions array +int npoints # Number of input points +bool debug # Print debugging information? + +int i, nlmax + +bool is_local_max() + +begin + nlmax = 0 + do i = 1, npoints { + if (is_local_max (i, data, npoints)) { + nlmax = nlmax + 1 + x[nlmax] = i + } + } + + if (debug) { + call printf (" Number of local minima/maxima found = %d.\n") + call pargi (nlmax) + } + + return (nlmax) +end + + +# IS_LOCAL_MAX -- Test a point to determine if it is a local minima/maximum. +# +# Indefinite points are ignored. + +bool procedure is_local_max (index, data, npoints) + +# Procedure parameters: +int index # Index to test for local minimum/maximum +real data[npoints] # Data values +int npoints # Number of points in the data vector + +int i, j, nright, nleft + +pointer y +int maxima +common /sort/ y, maxima + +begin + # INDEF points cannot be local minima/axima. + if (IS_INDEFR (data[index])) + return (FALSE) + + # Find the left and right indices where data values change and the + # number of points with the same value. Ignore INDEF points. + nleft = 0 + for (i = index - 1; i >= 1; i = i - 1) { + if (!IS_INDEFR (data[i])) { + if (data[i] != data[index]) + break + nleft = nleft + 1 + } + } + nright = 0 + for (j = index + 1; i <= npoints; j = j + 1) { + if (!IS_INDEFR (data[j])) { + if (data[j] != data[index]) + break + nright = nright + 1 + } + } + + # Test for failure to be a local minima/axima + if (maxima == YES) { + if ((i == 0) && (j == npoints)) { + return (FALSE) # Data is constant + } else if (i == 0) { + if (data[j] > data[index]) + return (FALSE) # Data increases to right + } else if (j == npoints) { + if (data[i] > data[index]) # Data increase to left + return (FALSE) + } else if ((data[i] > data[index]) || (data[j] > data[index])) { + return (FALSE) # Not a local maximum + } else if (!((nleft - nright == 0) || (nleft - nright == 1))) { + return (FALSE) # Not center of plateau + } + } else { + if ((i == 0) && (j == npoints)) { + return (FALSE) # Data is constant + } else if (i == 0) { + if (data[j] < data[index]) + return (FALSE) # Data decreases to right + } else if (j == npoints) { + if (data[i] < data[index]) # Data decrease to left + return (FALSE) + } else if ((data[i] < data[index]) || (data[j] < data[index])) { + return (FALSE) # Not a local maximum + } else if (!((nleft - nright == 0) || (nleft - nright == 1))) { + return (FALSE) # Not center of plateau + } + } + + # Point is a local minima/maxima + return (TRUE) +end + + +# FIND_THRESHOLD -- Find the peaks with positions satisfying threshold +# and contrast constraints. +# +# The input is the data array, data, and the peak positions array, x. +# The x array is resorted to the nthreshold peaks satisfying the constraints. +# The corresponding nthreshold data values are returned the y array. +# The number of peaks satisfying the constraints (nthreshold) is returned. + +int procedure find_threshold (data, x, y, npoints, contrast, threshold, debug) + +real data[ARB] # Input data values +real x[npoints] # Input/Output peak positions +real y[npoints] # Output peak data values +int npoints # Number of peaks input +real contrast # Contrast constraint +real threshold # Threshold constraint +bool debug # Print debugging information? + +int i, j, nthreshold +real minval, maxval, lcut + +pointer dummy +int maxima +common /sort/ dummy, maxima + +begin + # Set the y array to be the values at the peak positions. + do i = 1, npoints { + j = x[i] + y[i] = data[j] + } + + # Determine the min and max values of the peaks. + call alimr (y, npoints, minval, maxval) + + # Set the threshold based on the max of the absolute threshold and the + # contrast. Use arlt to set peaks below threshold to INDEF. + + if (maxima == YES) { + lcut = max (real (threshold), real (contrast * maxval)) + call arltr (y, npoints, lcut, INDEFR) + } else { + lcut = threshold + call argtr (y, npoints, lcut, INDEFR) + } + + if (debug) { + call printf (" Highest peak value = %g.\n") + call pargr (maxval) + call printf (" Peak cutoff threshold = %g.\n") + call pargr (lcut) + do i = 1, npoints { + if (IS_INDEFR (y[i])) { + j = x[i] + call printf ( + " Peak at column %d with value %g below threshold.\n") + call pargi (j) + call pargr (data[j]) + } + } + } + + # Determine the number of acceptable peaks & resort the x and y arrays. + nthreshold = 0 + do i = 1, npoints { + if (IS_INDEFR (y[i])) + next + nthreshold = nthreshold + 1 + x[nthreshold] = x[i] + y[nthreshold] = y[i] + } + + if (debug) { + call printf (" Number of peaks above the threshold = %d.\n") + call pargi (nthreshold) + } + + return (nthreshold) +end + +# FIND_ISOLATED -- Flag peaks which are within separation of a peak +# with a higher peak value. +# +# The peak positions, x, and their ranks, rank, are input. +# The rank array contains the indices of the peak positions in order from +# the highest peak value to the lowest peak value. Starting with +# highest rank (rank[1]) all peaks of lower rank within separation +# are marked by setting their positions to INDEF. The number of +# unflaged peaks is returned. + +int procedure find_isolated (x, rank, npoints, separation, debug) + +# Procedure parameters: +real x[npoints] # Positions of points +int rank[npoints] # Rank of peaks +int npoints # Number of peaks +int separation # Minimum allowed separation +bool debug # Print diagnostic information + +int i, j +int nisolated + +begin + # Eliminate close neighbors. The eliminated + # peaks are marked by setting their positions to INDEF. + nisolated = 0 + do i = 1, npoints { + if (IS_INDEFR (x[rank[i]])) + next + nisolated = nisolated + 1 + do j = i + 1, npoints { + if (IS_INDEFR (x[rank[j]])) + next + if (abs (x[rank[i]] - x[rank[j]]) < separation) { + if (debug) { + call printf ( + " Peak at column %d too near peak at column %d.\n") + call pargi (int (x[rank[j]])) + call pargi (int (x[rank[i]])) + } + x[rank[j]] = INDEFR + } + } + } + + if (debug) { + call printf (" Number of peaks separated by %d pixels = %d.\n") + call pargi (separation) + call pargi (nisolated) + } + + # Return number of isolated peaks. + return (nisolated) +end + + +# FIND_NMAX -- Select up to the nmax highest ranked peaks. +# +# The data values, data, peak positions, x, and their ranks, rank, are input. +# The data values are used only in printing debugging information. +# Peak positions previously eliminated are flaged by the value INDEF. +# The rank array contains the indices to the peak positions in order from +# the highest peak value to the lowest peak value. +# First all but the nmax highest ranked peaks (which have not been previously +# eliminated) are eliminated by marking their positions with the value INDEF. +# Then the remaining peaks are resorted to contain only the unflaged +# peaks and the number of such peaks is returned. + +int procedure find_nmax (data, x, rank, npoints, nmax, debug) + +real data[ARB] # Input data values +real x[npoints] # Peak positions +int rank[npoints] # Ranks of peaks +int npoints # Number of input peaks +int nmax # Max number of peaks to be selected +bool debug # Print debugging information? + +int i, j, npeaks + +begin + # Only mark peaks to reject if the number peaks is greater than nmax. + if (nmax < npoints) { + npeaks = 0 + do i = 1, npoints { + if (IS_INDEFR (x[rank[i]])) + next + npeaks = npeaks + 1 + if (npeaks > nmax) { + if (debug) { + j = x[rank[i]] + call printf ( + " Reject peak at column %d with rank %d and value %g.\n") + call pargi (j) + call pargi (i) + call pargr (data[j]) + } + x[rank[i]] = INDEFR + } + } + } + + # Eliminate INDEF points and determine the number of spectra found. + npeaks = 0 + do i = 1, npoints { + if (IS_INDEFR (x[i])) + next + npeaks = npeaks + 1 + x[npeaks] = x[i] + } + + return (npeaks) +end + + +# COMPARE -- Compare procedure for sort used in FIND_PEAKS. +# Larger values are indexed first. INDEF values are indexed last. + +int procedure compare (index1, index2) + +# Procedure parameters: +int index1 # Comparison index +int index2 # Comparison index + +pointer y +int maxima +common /sort/ y, maxima + +begin + # INDEF points are considered to be smallest/largest possible values. + if (maxima == YES) { + if (IS_INDEFR (Memr[y - 1 + index1])) + return (1) + else if (IS_INDEFR (Memr[y - 1 + index2])) + return (-1) + else if (Memr[y - 1 + index1] < Memr[y - 1 + index2]) + return (1) + else if (Memr[y - 1 + index1] > Memr[y - 1 + index2]) + return (-1) + else + return (0) + } else { + if (IS_INDEFR (Memr[y - 1 + index1])) + return (-1) + else if (IS_INDEFR (Memr[y - 1 + index2])) + return (1) + else if (Memr[y - 1 + index1] < Memr[y - 1 + index2]) + return (-1) + else if (Memr[y - 1 + index1] > Memr[y - 1 + index2]) + return (1) + else + return (0) + } +end diff --git a/noao/rv/rvidlines/reidentify.x b/noao/rv/rvidlines/reidentify.x new file mode 100644 index 00000000..ed8133b2 --- /dev/null +++ b/noao/rv/rvidlines/reidentify.x @@ -0,0 +1,609 @@ +include <error.h> +include <imhdr.h> +include <gset.h> +include "identify.h" + +define HELP "noao$onedspec/identify/identify.key" +define RVHELP "noao$rv/rvidlines/rvidlines.key" +define ICFITHELP "noao$lib/scr/idicgfit.key" +define PROMPT "identify options" + +define PAN 1 # Pan graph +define ZOOM 2 # Zoom graph + +# REIDENTIFY -- Reidentify features in an image. + +procedure reidentify (id) + +pointer id # ID pointer + +real wx, wy, wx2, wy2 +int wcs, key +char cmd[SZ_LINE] + +char newimage[SZ_FNAME] +int i, j, last, all, prfeature, nfeatures1, npeaks, ftype +double pix, fit, user, shift, pix_shift, z_shift, xg[10] +pointer peaks, label + +int clgcur(), scan(), nscan(), find_peaks(), errcode() +double id_center(), fit_to_pix(), id_fitpt() +double id_shift(), id_rms(), id_zshiftd(), id_zval() +errchk id_graph() + +define newim_ 10 +define newkey_ 20 +define beep_ 99 + +begin + # Initialize. + ID_GTYPE(id) = PAN + all = 0 + last = ID_CURRENT(id) + newimage[1] = EOS + ID_REFIT(id) = NO + ID_NEWFEATURES(id) = NO + ID_NEWCV(id) = NO + wy = INDEF + key = 'r' + + repeat { + prfeature = YES + if (all != 0) + all = mod (all + 1, 3) + + switch (key) { + case '?': # Print help + if (ID_TASK(id) == IDENTIFY) + call gpagefile (ID_GP(id), HELP, PROMPT) + else + call gpagefile (ID_GP(id), RVHELP, PROMPT) + case ':': # Process colon commands + if (cmd[1] == '/') + call gt_colon (cmd, ID_GP(id), ID_GT(id), ID_NEWGRAPH(id)) + else + call id_colon (id, cmd, newimage, prfeature) + case ' ': # Go to current feature + case '.': # Go to nearest feature + if (ID_NFEATURES(id) == 0) + goto beep_ + call id_nearest (id, double (wx)) + case '-': # Go to previous feature + if (ID_CURRENT(id) == 1) + goto beep_ + ID_CURRENT(id) = ID_CURRENT(id) - 1 + case '+', 'n': # Go to next feature + if (ID_CURRENT(id) == ID_NFEATURES(id)) + goto beep_ + ID_CURRENT(id) = ID_CURRENT(id) + 1 + case 'a': # Set all flag for next key + all = 1 + case 'b': # Mark blended features + i = 1 + fit = wx + xg[i] = fit_to_pix (id, fit) + call printf ("mark other components (exit with 'q'):") + while (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)!=EOF) { + if (key == 'q') + break + i = i + 1 + fit = wx + xg[i] = fit_to_pix (id, fit) + if (i == 10) + break + } + + call printf ("mark two background points: ") + j = clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) + j = clgcur ("cursor", wx2, wy2, wcs, key, cmd, SZ_LINE) + wx = fit_to_pix (id, double (wx)) + wx2 = fit_to_pix (id, double (wx2)) + + switch (ID_FTYPE(id)) { + case EMISSION: + ftype = GEMISSION + case ABSORPTION: + ftype = GABSORPTION + default: + ftype = ID_FTYPE(id) + } + iferr (call id_gcenter (id, xg, i, wx, wy, wx2, wy2, + ID_FWIDTH(id), ftype, YES)) { + call erract (EA_WARN) + prfeature = NO + goto newkey_ + } + + do j = 1, i { + pix = xg[j] + fit = id_fitpt (id, pix) + user = fit + call id_newfeature (id, pix, fit, user, 1.0D0, + ID_FWIDTH(id), ftype, NULL) + USER(id,ID_CURRENT(id)) = INDEFD + call id_match (id, FIT(id,ID_CURRENT(id)), + USER(id,ID_CURRENT(id)), + Memi[ID_LABEL(id)+ID_CURRENT(id)-1], + ID_MATCH(id)) + call id_mark (id, ID_CURRENT(id)) + call printf ("%10.2f %10.8g ") + call pargd (PIX(id,ID_CURRENT(id))) + call pargd (FIT(id,ID_CURRENT(id))) + if (ID_REDSHIFT(id) != 0.) { + call printf ("%10.8g ") + call pargd ( + id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0)) + } + call printf ("(%10.8g %s): ") + call pargd (USER(id,ID_CURRENT(id))) + label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + call gargwrd (cmd, SZ_LINE) + i = nscan() + if (i > 0) { + USER(id,ID_CURRENT(id)) = user + fit = id_zshiftd (id, user, 1) + call id_match (id, fit, USER(id,ID_CURRENT(id)), + Memi[ID_LABEL(id)+ID_CURRENT(id)-1], + ID_MATCH(id)) + } + if (i > 1) { + call reset_scan () + call gargd (user) + call gargstr (cmd, SZ_LINE) + call id_label (cmd, + Memi[ID_LABEL(id)+ID_CURRENT(id)-1]) + } + } + } + case 'c': # Recenter features + if (all != 0) { + for (i = 1; i <= ID_NFEATURES(id); i = i + 1) { + call gseti (ID_GP(id), G_PLTYPE, 0) + call id_mark (id, i) + call gseti (ID_GP(id), G_PLTYPE, 1) + FWIDTH(id,i) = ID_FWIDTH(id) + PIX(id,i) = id_center (id, PIX(id,i), 1, FWIDTH(id,i), + FTYPE(id,i), NO) + if (!IS_INDEFD (PIX(id,i))) { + FIT(id,i) = id_fitpt (id, PIX(id,i)) + call id_mark (id, i) + } else { + call id_delete (id, i) + i = i - 1 + } + } + ID_NEWFEATURES(id) = YES + } else { + if (ID_NFEATURES(id) < 1) + goto beep_ + call id_nearest (id, double (wx)) + pix = PIX(id,ID_CURRENT(id)) + pix = id_center (id, pix, 1, ID_FWIDTH(id), + FTYPE(id,ID_CURRENT(id)), NO) + if (!IS_INDEFD (pix)) { + call gseti (ID_GP(id), G_PLTYPE, 0) + call id_mark (id, ID_CURRENT(id)) + PIX(id,ID_CURRENT(id)) = pix + FWIDTH(id,ID_CURRENT(id)) = ID_FWIDTH(id) + FIT(id,ID_CURRENT(id)) = id_fitpt (id, pix) + call gseti (ID_GP(id), G_PLTYPE, 1) + call id_mark (id, ID_CURRENT(id)) + ID_NEWFEATURES(id) = YES + } else { + call printf ("Centering failed\n") + prfeature = NO + } + } + case 'd': # Delete features + if (all != 0) { + ID_NFEATURES(id) = 0 + ID_CURRENT(id) = 0 + ID_NEWFEATURES(id) = YES + ID_NEWGRAPH(id) = YES + } else { + if (ID_NFEATURES(id) < 1) + goto beep_ + call id_nearest (id, double (wx)) + call gseti (ID_GP(id), G_PLTYPE, 0) + call id_mark (id, ID_CURRENT(id)) + call gseti (ID_GP(id), G_PLTYPE, 1) + call id_delete (id, ID_CURRENT(id)) + ID_CURRENT(id) = min (ID_NFEATURES(id), ID_CURRENT(id)) + last = 0 + } + case 'f': # Fit dispersion function + if (ID_TASK(id) == IDENTIFY) { + call id_dofit (id, YES) + } else { + call id_velocity (id, YES) + prfeature = NO + } + case 'g': # Fit shift + if (ID_TASK(id) == RVIDLINES) + goto beep_ + + call id_doshift (id, YES) + prfeature = NO + case 'i': # Initialize + call dcvfree (ID_CV(id)) + ID_SHIFT(id) = 0. + ID_REDSHIFT(id) = 0. + ID_NEWCV(id) = YES + ID_NFEATURES(id) = 0 + ID_CURRENT(id) = 0 + ID_NEWFEATURES(id) = YES + ID_NEWGRAPH(id) = YES + case 'j', 'k', 'o': + call printf ("Command not available in REIDENTIFY") + prfeature = NO + case 'l': # Find features from line list + if (ID_TASK(id) == IDENTIFY) { + if (ID_NFEATURES(id) >= 2) + call id_dofit (id, NO) + if (ID_NEWCV(id) == YES) { + iferr (call id_fitdata(id)) + ; + call id_fitfeatures(id) + ID_NEWCV(id) = NO + } + call id_linelist (id) + if (ID_NEWFEATURES(id) == YES) + ID_REFIT(id) = YES + } else { + call id_velocity (id, NO) + call id_linelist(id) + if (ID_NEWFEATURES(id) == YES) { + call id_velocity (id, NO) + ID_NEWGRAPH(id) = YES + } + } + case 'm': # Mark new feature + fit = wx + pix = fit_to_pix (id, fit) + pix = id_center (id, pix, 1, ID_FWIDTH(id), ID_FTYPE(id), YES) + if (IS_INDEFD (pix)) + goto beep_ + fit = id_fitpt (id, pix) + user = fit + call id_newfeature (id, pix, fit, user, 1.0D0, ID_FWIDTH(id), + ID_FTYPE(id), NULL) + USER(id,ID_CURRENT(id)) = INDEFD + call id_match (id, FIT(id,ID_CURRENT(id)), + USER(id,ID_CURRENT(id)), + Memi[ID_LABEL(id)+ID_CURRENT(id)-1], + ID_MATCH(id)) + call id_mark (id, ID_CURRENT(id)) + call printf ("%10.2f %10.8g ") + call pargd (PIX(id,ID_CURRENT(id))) + call pargd (FIT(id,ID_CURRENT(id))) + if (ID_REDSHIFT(id) != 0.) { + call printf ("[%10.8g] ") + call pargd ( + id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0)) + } + call printf ("(%10.8g %s): ") + call pargd (USER(id,ID_CURRENT(id))) + label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + call gargwrd (cmd, SZ_LINE) + i = nscan() + if (i > 0) { + USER(id,ID_CURRENT(id)) = user + call id_match (id, user, USER(id,ID_CURRENT(id)), + Memi[ID_LABEL(id)+ID_CURRENT(id)-1], + ID_MATCH(id)) + } + if (i > 1) { + call reset_scan () + call gargd (user) + call gargstr (cmd, SZ_LINE) + call id_label (cmd, Memi[ID_LABEL(id)+ID_CURRENT(id)-1]) + } + } + case 'p': # Switch to pan mode + if (ID_GTYPE(id) != PAN) { + ID_GTYPE(id) = PAN + ID_NEWGRAPH(id) = YES + } + case 'q': # Exit loop + break + case 'r': # Redraw the graph + ID_NEWGRAPH(id) = YES + case 's', 'x': # Shift or correlate features + if (ID_TASK(id) == RVIDLINES) + goto beep_ + + # Get coordinate shift. + switch (key) { + case 's': + call printf ("User coordinate (%10.8g): ") + call pargr (wx) + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + if (nscan() == 1) + shift = wx - user + } else + shift = 0. + case 'x': + if (ID_NFEATURES(id) > 5) + shift = id_shift (id) + else + goto beep_ + } + + ID_NEWFEATURES(id) = YES + ID_NEWCV(id) = YES + ID_NEWGRAPH(id) = YES + prfeature = NO + + if (ID_NFEATURES(id) < 1) { + call printf ("User coordinate shift=%5f\n") + call pargd (shift) + ID_SHIFT(id) = ID_SHIFT(id) + shift + goto newkey_ + } + + # Recenter features. + pix_shift = 0. + z_shift = 0. + nfeatures1 = ID_NFEATURES(id) + + j = 0. + do i = 1, ID_NFEATURES(id) { + pix = fit_to_pix (id, FIT(id,i) + shift) + pix = id_center (id, pix, 1, FWIDTH(id,i), FTYPE(id,i), NO) + if (IS_INDEFD (pix)) { + if (ID_CURRENT(id) == i) + ID_CURRENT(id) = i + 1 + next + } + fit = id_fitpt (id, pix) + + pix_shift = pix_shift + pix - PIX(id,i) + if (FIT(id,i) != 0.) + z_shift = z_shift + id_zval (id, fit, FIT(id,i)) + + j = j + 1 + PIX(id,j) = pix + FIT(id,j) = FIT(id,i) + USER(id,j) = USER(id,i) + WTS(id,j) = WTS(id,i) + FWIDTH(id,j) = FWIDTH(id,i) + FTYPE(id,j) = FTYPE(id,i) + if (ID_CURRENT(id) == i) + ID_CURRENT(id) = j + } + if (j != ID_NFEATURES(id)) { + ID_NFEATURES(id) = j + ID_CURRENT(id) = min (ID_CURRENT(id), ID_NFEATURES(id)) + } + + if (ID_NFEATURES(id) < 1) { + call printf ("User coordinate shift=%5f") + call pargd (shift) + call printf (", No features found during recentering\n") + ID_SHIFT(id) = ID_SHIFT(id) + shift + goto newkey_ + } + + # Adjust shift. + pix = ID_SHIFT(id) + call id_doshift (id, NO) + call id_fitfeatures (id) + + # Print results. + call printf ("Recentered=%d/%d") + call pargi (ID_NFEATURES(id)) + call pargi (nfeatures1) + call printf ( + ", pixel shift=%.2f, user shift=%5f, z=%7.3g, rms=%5g\n") + call pargd (pix_shift / ID_NFEATURES(id)) + call pargd (pix - ID_SHIFT(id)) + call pargd (z_shift / ID_NFEATURES(id)) + call pargd (id_rms(id)) + case 't': # Move the current feature + if (ID_CURRENT(id) < 1) + goto beep_ + pix = fit_to_pix (id, double (wx)) + call gseti (ID_GP(id), G_PLTYPE, 0) + call id_mark (id, ID_CURRENT(id)) + PIX(id,ID_CURRENT(id)) = pix + FIT(id,ID_CURRENT(id)) = id_fitpt (id, pix) + call gseti (ID_GP(id), G_PLTYPE, 1) + call id_mark (id, ID_CURRENT(id)) + ID_NEWFEATURES(id) = YES + case 'u': # Set user coordinate + if (ID_NFEATURES(id) < 1) + goto beep_ + call printf ("%10.2f %10.8g ") + call pargd (PIX(id,ID_CURRENT(id))) + call pargd (FIT(id,ID_CURRENT(id))) + if (ID_REDSHIFT(id) != 0.) { + call printf ("[%10.8g] ") + call pargd (id_zshiftd (id, FIT(id,ID_CURRENT(id)), 0)) + } + call printf ("(%10.8g %s): ") + call pargd (USER(id,ID_CURRENT(id))) + label = Memi[ID_LABEL(id)+ID_CURRENT(id)-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + call gargwrd (cmd, SZ_LINE) + i = nscan() + if (i > 0) { + USER(id,ID_CURRENT(id)) = user + ID_NEWFEATURES(id) = YES + } + if (i > 1) { + call reset_scan () + call gargd (user) + call gargstr (cmd, SZ_LINE) + call id_label (cmd, Memi[ID_LABEL(id)+ID_CURRENT(id)-1]) + } + } + call printf ("Weight (%g): ") + call pargd (WTS(id,ID_CURRENT(id))) + call flush (STDOUT) + if (scan() != EOF) { + call gargd (user) + if (nscan() > 0) + WTS(id,ID_CURRENT(id)) = user + } + case 'w': # Window graph + call gt_window (ID_GT(id), ID_GP(id), "cursor", ID_NEWGRAPH(id)) + case 'y': # Find peaks + call malloc (peaks, ID_NPTS(id), TY_REAL) + npeaks = find_peaks (IMDATA(id,1), Memr[peaks], ID_NPTS(id), 0., + int (ID_MINSEP(id)), 0, ID_MAXFEATURES(id), 0., false) + for (j = 1; j <= ID_NFEATURES(id); j = j + 1) { + for (i = 1; i <= npeaks; i = i + 1) { + if (!IS_INDEF (Memr[peaks+i-1])) { + pix = Memr[peaks+i-1] + if (abs (pix - PIX(id,j)) < ID_MINSEP(id)) + Memr[peaks+i-1] = INDEF + } + } + } + for (i = 1; i <= npeaks; i = i + 1) { + if (IS_INDEF(Memr[peaks+i-1])) + next + pix = Memr[peaks+i-1] + pix = id_center(id, pix, 1, ID_FWIDTH(id), ID_FTYPE(id), NO) + if (IS_INDEFD (pix)) + next + fit = id_fitpt (id, pix) + user = INDEFD + call id_match (id, fit, user, label, ID_MATCH(id)) + call id_newfeature (id, pix, fit, user, 1.0D0, + ID_FWIDTH(id), ID_FTYPE(id), label) + call id_mark (id, ID_CURRENT(id)) + } + call mfree (peaks, TY_REAL) + case 'z': # Go to zoom mode + if (ID_NFEATURES(id) < 1) + goto beep_ + if (ID_GTYPE(id) == PAN) + ID_NEWGRAPH(id) = YES + ID_GTYPE(id) = ZOOM + call id_nearest (id, double (wx)) + case 'I': + call fatal (0, "Interrupt") + default: +beep_ call printf ("\007") + } + +newkey_ + # Set update flag if anything has changed. + if ((ID_NEWFEATURES(id) == YES) || (ID_NEWCV(id) == YES)) + ID_NEWDBENTRY(id) = YES + + # If a new image exit loop, update database, and start over. + if (newimage[1] != EOS) { + call printf ("Can't change image in REIDENTIFY") + newimage[1] = EOS + prfeature = NO + } + + # Refit dispersion function + if (ID_REFIT(id) == YES) { + call id_dofit (id, NO) + ID_REFIT(id) = NO + } + + # If there is a new dispersion solution evaluate the coordinates + if (ID_NEWCV(id) == YES) { + iferr (call id_fitdata (id)) + ; + call id_fitfeatures (id) + ID_NEWCV(id) = NO + } + + # Draw new graph in zoom mode if current feature has changed. + if ((ID_GTYPE(id) == ZOOM) && (last != ID_CURRENT(id))) + ID_NEWGRAPH(id) = YES + + # Draw new graph. + if (ID_NEWGRAPH(id) == YES) { + call id_graph (id, ID_GTYPE(id)) + ID_NEWGRAPH(id) = NO + } + + # Set cursor and print status of current feature (unless canceled). + if (ID_CURRENT(id) > 0) { + if (IS_INDEF (wy)) { + i = max (1, min (ID_NPTS(id), int (PIX(id,ID_CURRENT(id))))) + wy = IMDATA(id,i) + } + + call gscur (ID_GP(id), real (FIT(id,ID_CURRENT(id))), wy) + if (errcode() == OK && prfeature == YES) { + i = ID_CURRENT(id) + pix = PIX(id,i) + fit = FIT(id,i) + user = USER(id,i) + if (ID_TASK(id) == IDENTIFY) { + if (IS_INDEFD(user)) + shift = INDEF + else + shift = fit - user + call printf ("%10.2f %10.8g %10.8g %10.3g %s") + call pargd (pix) + call pargd (fit) + call pargd (user) + call pargd (shift) + label = Memi[ID_LABEL(id)+i-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + } else { + if (IS_INDEFD(user)) + shift = INDEF + else { + shift = id_zval (id, fit, user) - ID_REDSHIFT(id) + if (abs (shift) < 0.01) + shift = shift * VLIGHT + } + call printf ("%10.2f %10.8g %10.8g %10.8g %10.4g %s") + call pargd (pix) + call pargd (fit) + call pargd (id_zshiftd (id, fit, 0)) + call pargd (user) + if (IS_INDEFD(user)) + call pargd (INDEFD) + else + call pargd (shift) + label = Memi[ID_LABEL(id)+i-1] + if (label != NULL) + call pargstr (Memc[label]) + else + call pargstr ("") + } + } + } + + # Print delayed error message + if (errcode() != OK) + call erract (EA_WARN) + + last = ID_CURRENT(id) + } until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF) +end diff --git a/noao/rv/rvidlines/rvidlines.key b/noao/rv/rvidlines/rvidlines.key new file mode 100644 index 00000000..cfccacde --- /dev/null +++ b/noao/rv/rvidlines/rvidlines.key @@ -0,0 +1,100 @@ + RVIDLINES HELP + + +STATUS LINE + +The status line gives the pixel position, observed wavelength, +rest wavelength, user wavelength, velocity residual, and optional +line identification: + + pixel observed rest user Vresidual [identification] + + +CURSOR KEY SUMMARY + +? Help l Match list w Window graph +a Affect all features m Mark feature y Find peaks +b Deblend n Next feature z Zoom graph +c Center feature(s) o Go to line + Next feature +d Delete feature(s) p Pan graph - Previous feature +f Fit velocity q Quit . Nearest feature +i Initialize r Redraw graph I Interrupt +j Preceding line t Reset position +k Next line u Enter coordinate + + +COLON COMMAND SUMMARY + +:show [file] :features [file] :coordlist [file] +:cradius [value] :threshold [value] :database [file] +:ftype [type] :fwidth [value] :image [image] +:labels [type] :match [value] :maxfeatures [value] +:minsep [value] :read [image [ap]] :write [image [ap]] +:zwidth [value] :add [image [ap]] + + +CURSOR KEYS + +? Clear the screen and print menu of options +a Apply next (c)enter or (d)elete operation to (a)ll features +b Mark and de(b)lend features by Gaussian fitting +c (C)enter the feature nearest the cursor +d (D)elete the feature nearest the cursor +f (F)it a redshift and velocity from the fitted and user coordinates +i (I)nitialize (delete features and coordinate fit) +j Go to the preceding image line or column in a 2D or multispec image +k Go to the next image line or column in a 2D or multispec image +l Match coordinates in the coordinate (l)ist +m (M)ark a new feature near the cursor and enter coordinate and label +n Move the cursor or zoom to the (n)ext feature (same as +) +o Go to the specified image line or column in a 2D or multispec image +p (P)an to user defined window after (z)ooming on a feature +q (Q)uit and continue with next image (also carriage return) +r (R)edraw the graph +t Reset the position of a feature without centering +u Enter a new (u)ser coordinate and label for the current feature +w (W)indow the graph. Use '?' to window prompt for more help. +y Automatically find "maxfeatures" strongest peaks and identify them +z (Z)oom on the feature nearest the cursor +. Move the cursor or zoom to the feature nearest the cursor ++ Move the cursor or zoom to the next feature +- Move the cursor or zoom to the previous feature +I Interrupt task and exit immediately. Database information is not saved. + + +COLON COMMANDS + +The parameters are listed or set with the following commands which may be +abbreviated. To list the value of a parameter type the command alone. + +:show file Show the values of all the parameters +:features file Write feature list to file (default is STDOUT) + +:coordlist file Coordinate list file +:cradius value Centering radius in pixels +:threshold value Detection threshold for feature centering +:database name Database for recording feature records +:ftype value Feature type + (emission|absorption|gemission|gabsorption) +:fwidth value Feature width in pixels +:image imagename Set a new image or show the current image +:labels value Feature label type + (none|index|pixel|coords|user|both) +:match value Coordinate list matching distance +:maxfeatures value Maximum number of features automatically found +:minsep value Minimum separation allowed between features +:read name ap Read a record from the database + (name and ap default to the current spectrum) +:write name ap Write a record to the database + (name and ap default to the current spectrum) +:add name ap Add features from the database + (name and ap default to the current spectrum) +:zwidth value Zoom width in user units + +Labels: + none - No labels + index - Sequential numbers in order of increasing pixel position + pixel - Pixel coordinates + coords - User coordinates such as wavelength + user - User labels + both - Combination of coords and user diff --git a/noao/rv/rvidlines/t_identify.x b/noao/rv/rvidlines/t_identify.x new file mode 100644 index 00000000..aac9596f --- /dev/null +++ b/noao/rv/rvidlines/t_identify.x @@ -0,0 +1,108 @@ +include <mach.h> +include <pkg/gtools.h> +include "identify.h" + +# T_IDENTIFY -- Identify features and determine dispersion solutions + +procedure t_identify () + +begin + call identify (IDENTIFY) +end + + +# T_RVIDLINES -- Identify features and determine radial velocities + +procedure t_rvidlines () + +begin + call identify (RVIDLINES) +end + + +# IDENTIFY -- Initialize task + +procedure identify (taskid) + +int taskid #I Task ID + +int list, clscan(), clgeti(), clgwrd(), nscan(), imtopenp(), imtgetim() +real clgetr() +pointer sp, str, id, gt_init() + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # Allocate the basic data structure. + call id_init (id, taskid) + + # Get task parameters. + list = imtopenp ("images") + if (clscan ("nsum") != EOF) { + call gargi (ID_NSUM(id,1)) + call gargi (ID_NSUM(id,2)) + if (nscan() == 0) + call error (1, "Error in 'nsum' parameter") + if (nscan() == 1) + ID_NSUM(id,2) = ID_NSUM(id,1) + ID_NSUM(id,1) = max (1, ID_NSUM(id,1)) + ID_NSUM(id,2) = max (1, ID_NSUM(id,2)) + } + ID_MAXFEATURES(id) = clgeti ("maxfeatures") + ID_MINSEP(id) = clgetr ("minsep") + ID_MATCH(id) = clgetr ("match") + ID_ZWIDTH(id) = clgetr ("zwidth") + ID_FTYPE(id) = clgwrd ("ftype", Memc[str], SZ_LINE, FTYPES) + ID_FWIDTH(id) = clgetr ("fwidth") + ID_CRADIUS(id) = clgetr ("cradius") + ID_THRESHOLD(id) = clgetr ("threshold") + call clgstr ("database", Memc[ID_DATABASE(id)], SZ_FNAME) + call clgstr ("coordlist", Memc[ID_COORDLIST(id)], SZ_FNAME) + ID_LABELS(id) = 1 + + # Initialize features data structure. + ID_GT(id) = gt_init() + call gt_sets (ID_GT(id), GTTYPE, "line") + ID_CV(id) = NULL + ID_CURRENT(id) = 0 + ID_SHIFT(id) = 0. + ID_REDSHIFT(id) = 0. + + # Initialize ICFIT + call ic_open (ID_IC(id)) + if (ID_TASK(id) == IDENTIFY) { + call clgstr ("function", Memc[str], SZ_LINE) + call ic_pstr (ID_IC(id), "function", Memc[str]) + call ic_puti (ID_IC(id), "order", clgeti ("order")) + call clgstr ("sample", Memc[str], SZ_LINE) + call ic_pstr (ID_IC(id), "sample", Memc[str]) + call ic_puti (ID_IC(id), "naverage", 1) + call ic_puti (ID_IC(id), "niterate", clgeti ("niterate")) + call ic_putr (ID_IC(id), "low", clgetr ("low_reject")) + call ic_putr (ID_IC(id), "high", clgetr ("high_reject")) + call ic_putr (ID_IC(id), "grow", clgetr ("grow")) + call ic_pstr (ID_IC(id), "xlabel", "Feature positions") + call ic_pstr (ID_IC(id), "xunits", "pixels") + call ic_pstr (ID_IC(id), "ylabel", "") + call ic_pkey (ID_IC(id), 1, 'y', 'x') + call ic_pkey (ID_IC(id), 2, 'y', 'v') + call ic_pkey (ID_IC(id), 3, 'y', 'r') + call ic_pkey (ID_IC(id), 4, 'y', 'd') + call ic_pkey (ID_IC(id), 5, 'y', 'n') + call ic_puti (ID_IC(id), "key", 3) + } + + # Get the line list. + call id_mapll (id) + + # Expand the image template and identify features in each image. + while (imtgetim (list, Memc[ID_IMAGE(id)], SZ_FNAME) != EOF) + call id_identify (id) + + # Finish up. + call smw_daxis (NULL, NULL, 0, 0, 0) + call id_free (id) + call imtclose (list) + call sfree (sp) +end diff --git a/noao/rv/rvidlines/t_reidentify.x b/noao/rv/rvidlines/t_reidentify.x new file mode 100644 index 00000000..5a724325 --- /dev/null +++ b/noao/rv/rvidlines/t_reidentify.x @@ -0,0 +1,1092 @@ +include <error.h> +include <fset.h> +include <gset.h> +include <pkg/gtools.h> +include <smw.h> +include "identify.h" + +define ICFITHELP "noao$lib/scr/idicgfit.key" +define VLIGHT 2.997925e5 # Speed of light, Km/sec + + +# T_REIDENTIFY -- Reidentify features starting from reference features. + +procedure t_reidentify () + +begin + call reiden (IDENTIFY) +end + + +# T_RVREIDLINES -- Reidentify lines for radial velocities + +procedure t_rvreidlines () + +begin + call reiden (RVIDLINES) +end + + +# REIDEN -- Reidentify features starting from reference features. +# A reference spectrum is specified and the same features are identified +# in other images. Some lines may be lost due to bad centering. Additional +# lines may be excluded from a new fit to the dispersion function. Instead +# of refitting the dispersion function the user may elect to determine only +# a shift in the reference dispersion function. Additional features may +# be added given a coordinate list. +# +# In 2D images a starting line or column is selected. A number of lines +# or columns may be averaged before identifying features. If a positive step +# size is given then additional lines or columns may be reidentified in +# the reference image. This may be done either by tracing or by reidentifying +# starting from the same reference features. Reidentification between images +# is done by taking the same line or column from the reference image. +# +# Multispec format images are matched by aperture number and the spectra +# need not be in the same order in each image. + +procedure reiden (taskid) + +int taskid #I Task ID + +pointer reference # Reference image +int list # List of images +char ans[3] # Interactive? + +int i, fd, nlogfd +pointer sp, logfile, str, id, logfd, pd + +int clscan(), clgeti(), clpopnu(), clgfil(), clgwrd(), btoi() +int nscan(), open(), nowhite(), imtopenp(), imtgetim() +bool clgetb() +real clgetr() +pointer gopen(), gt_init() + +begin + call smark (sp) + call salloc (reference, SZ_FNAME, TY_CHAR) + call salloc (logfile, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Allocate the basic data structures. + call id_init (id, taskid) + call ic_open (ID_IC(id)) + + # Get task parameters. + call clgstr ("reference", Memc[reference], SZ_FNAME) + list = imtopenp ("images") + i = nowhite (Memc[reference], Memc[reference], SZ_FNAME) + + if (ID_TASK(id) == IDENTIFY) + ID_REFIT(id) = btoi (clgetb ("refit")) + else + ID_REFIT(id) = 3 + + if (clscan ("nsum") != EOF) { + call gargi (ID_NSUM(id,1)) + call gargi (ID_NSUM(id,2)) + if (nscan() == 0) + call error (1, "Error in 'nsum' parameter") + if (nscan() == 1) + ID_NSUM(id,2) = ID_NSUM(id,1) + ID_NSUM(id,1) = max (1, ID_NSUM(id,1)) + ID_NSUM(id,2) = max (1, ID_NSUM(id,2)) + } + ID_MAXFEATURES(id) = clgeti ("maxfeatures") + ID_MINSEP(id) = clgetr ("minsep") + ID_MATCH(id) = clgetr ("match") + if (ID_TASK(id) == IDENTIFY) { + ID_ZWIDTH(id) = clgetr ("identify.zwidth") + ID_FTYPE(id) = clgwrd ("identify.ftype", Memc[str], SZ_LINE, FTYPES) + ID_FWIDTH(id) = clgetr ("identify.fwidth") + } else { + ID_ZWIDTH(id) = clgetr ("rvidlines.zwidth") + ID_FTYPE(id) = clgwrd ("rvidlines.ftype", Memc[str], SZ_LINE,FTYPES) + ID_FWIDTH(id) = clgetr ("rvidlines.fwidth") + } + ID_CRADIUS(id) = clgetr ("cradius") + ID_THRESHOLD(id) = clgetr ("threshold") + call clgstr ("database", Memc[ID_DATABASE(id)], SZ_FNAME) + call clgstr ("coordlist", Memc[ID_COORDLIST(id)], SZ_FNAME) + ID_LABELS(id) = 1 + + call id_mapll (id) + ID_LOGFILES(id) = clpopnu ("logfiles") + + switch (clgwrd ("interactive", ans, SZ_FNAME, "|no|yes|NO|YES|")) { + case 1, 3: + call strcpy ("NO", ans, 3) + ID_GP(id) = NULL + case 2, 4: + # Open graphics + call clgstr ("graphics", Memc[logfile], SZ_FNAME) + ID_GP(id) = gopen (Memc[logfile], NEW_FILE+AW_DEFER, STDGRAPH) + if (ID_TASK(id) == IDENTIFY) { + call ic_pstr (ID_IC(id), "help", ICFITHELP) + call ic_pstr (ID_IC(id), "xlabel", "Feature positions") + call ic_pstr (ID_IC(id), "xunits", "pixels") + call ic_pstr (ID_IC(id), "ylabel", "") + call ic_pkey (ID_IC(id), 1, 'y', 'x') + call ic_pkey (ID_IC(id), 2, 'y', 'v') + call ic_pkey (ID_IC(id), 3, 'y', 'r') + call ic_pkey (ID_IC(id), 4, 'y', 'd') + call ic_pkey (ID_IC(id), 5, 'y', 'n') + call ic_puti (ID_IC(id), "key", 3) + } + } + + # Open log and plot files. + nlogfd = 0 + if (clgetb ("verbose")) { + nlogfd = 1 + call malloc (logfd, nlogfd, TY_INT) + Memi[logfd] = STDOUT + } + while (clgfil (ID_LOGFILES(id), Memc[logfile], SZ_FNAME) != EOF) { + fd = open (Memc[logfile], APPEND, TEXT_FILE) + call fseti (fd, F_FLUSHNL, YES) + nlogfd = nlogfd + 1 + if (nlogfd == 1) + call malloc (logfd, nlogfd, TY_INT) + else + call realloc (logfd, nlogfd, TY_INT) + Memi[logfd+nlogfd-1] = fd + } + call ri_loghdr (id, Memc[reference], Memi[logfd], nlogfd, 1) + + pd = NULL + if (ID_TASK(id) == IDENTIFY) { + call clgstr ("plotfile", Memc[logfile], SZ_FNAME) + if (nowhite (Memc[logfile], Memc[logfile], SZ_FNAME) > 0) { + fd = open (Memc[logfile], APPEND, BINARY_FILE) + pd = gopen ("stdvdm", NEW_FILE, fd) + } + } + + ID_GT(id) = gt_init() + call gt_sets (ID_GT(id), GTTYPE, "line") + + # Get and trace the reference solutions. + call ri_reference (id, Memc[reference], ans, Memi[logfd], nlogfd, pd) + + # Expand the image template and reidentify features. + while (imtgetim (list, Memc[ID_IMAGE(id)], SZ_FNAME) != EOF) + call ri_image (id, Memc[reference], Memc[ID_IMAGE(id)], ans, + Memi[logfd], nlogfd, pd) + + # Finish up. + if (nlogfd > 0) { + do i = 1, nlogfd + call close (Memi[logfd+i-1]) + call mfree (logfd, TY_INT) + } + + if (ID_GP(id) != NULL) + call gclose (ID_GP(id)) + if (pd != NULL) { + call gclose (pd) + call close (fd) + } + call clpcls (ID_LOGFILES(id)) + call imtclose (list) + call id_free (id) + call smw_daxis (NULL, NULL, 0, 0, 0) + call sfree (sp) +end + + +# RI_REFERENCE -- Set reference features. Trace if needed. + +procedure ri_reference (id, reference, ans, logfd, nlogfd, pd) + +pointer id # ID pointer +char reference[ARB] # Reference image +char ans[3] # Interactive? +int logfd[ARB] # Logfiles +int nlogfd # Number of logfiles +pointer pd # Plot file pointer + +int step[2] +double shift[2] +int nreid +bool override +bool trace + +int i, start[2], line[2], loghdr +double fit_shift[2] +pointer ic, ic1 +bool clgetb() +int clscan(), clgeti(), nscan(), id_gid(), id_dbcheck() +errchk id_dbread + +begin + # Open the image and return if there is an error. + call strcpy (reference, Memc[ID_IMAGE(id)], SZ_FNAME) + iferr (call id_map (id)) { + call erract (EA_WARN) + return + } + + # Set parameters + start[1] = ID_LINE(id,1) + start[2] = ID_LINE(id,2) + if (clscan ("step") != EOF) { + call gargi (step[1]) + call gargi (step[2]) + if (nscan() == 0) + call error (1, "Error in 'step' parameter") + if (nscan() == 1) + step[2] = step[1] + } + if (clscan ("shift") != EOF) { + call gargd (shift[1]) + call gargd (shift[2]) + if (nscan() == 0) + call error (1, "Error in 'shift' parameter") + if (nscan() == 1) + shift[2] = shift[1] + } + nreid = max (1, ID_NFEATURES(id) - clgeti ("nlost")) + override = clgetb ("override") + trace = clgetb ("trace") + + # Get and save the reference database entry. + call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO, NO) + call id_saveid (id, ID_LINE(id,1)) + + # Get and save other entries. + if (!override) { + for (line[2]=start[2]; line[2]>0; line[2]=line[2]-step[2]) { + ID_LINE(id,2) = line[2] + ID_AP(id,2) = line[2] + for (line[1]=start[1]; line[1]>0; line[1]=line[1]-step[1]) { + if (line[1]==start[1] && line[2]==start[2]) + next + ID_LINE(id,1) = line[1] + ID_AP(id,1) = line[1] + if (ID_APS(id) != NULL) + ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1] + ifnoerr ( + call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), + NO, NO)) { + call id_saveid (id, ID_LINE(id,1)) + } + } + for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1); + line[1]=line[1]+step[1]) { + ID_LINE(id,1) = line[1] + ID_AP(id,1) = line[1] + if (ID_APS(id) != NULL) + ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1] + ifnoerr (call id_dbread (id, Memc[ID_IMAGE(id)], + ID_AP(id,1), NO, NO)) { + call id_saveid (id, ID_LINE(id,1)) + } + } + } + for (line[2]=start[2]+step[2]; line[2]<=ID_MAXLINE(id,2); + line[2]=line[2]+step[2]) { + ID_LINE(id,2) = line[2] + ID_AP(id,2) = line[2] + for (line[1]=start[1]-step[1]; line[1]>0; + line[1]=line[1]-step[1]) { + ID_LINE(id,1) = line[1] + ID_AP(id,1) = line[1] + if (ID_APS(id) != NULL) + ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1] + ifnoerr ( + call id_dbread (id, Memc[ID_IMAGE(id)], ID_AP(id,1), + NO, NO)) { + call id_saveid (id, ID_LINE(id,1)) + } + } + for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1); + line[1]=line[1]+step[1]) { + ID_LINE(id,1) = line[1] + ID_AP(id,1) = line[1] + if (ID_APS(id) != NULL) + ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1] + ifnoerr (call id_dbread (id, Memc[ID_IMAGE(id)], + ID_AP(id,1), NO, NO)) { + call id_saveid (id, ID_LINE(id,1)) + } + } + } + } + + # Reidentify. + loghdr = 2 + ic = ID_IC(id) + if (ans[1] == 'N') + ic1 = ic + else { + call ic_open (ic1) + call ic_copy (ic, ic1) + } + + fit_shift[2] = shift[2] + for (line[2]=start[2]; line[2]>0; line[2]=line[2]-step[2]) { + ID_LINE(id,2) = line[2] + ID_AP(id,2) = line[2] + ID_IC(id) = ic + + if (IS_INDEFD(shift[2])) + fit_shift[2] = INDEFD + else { + if (!trace) + fit_shift[2] = fit_shift[2] - shift[2] + else + fit_shift[2] = -shift[2] + } + + fit_shift[1] = fit_shift[2] + for (line[1]=start[1]; line[1]>0; line[1]=line[1]-step[1]) { + if (line[1]==start[1] && line[2]==start[2]) + next + ID_LINE(id,1) = line[1] + ID_AP(id,1) = line[1] + ID_IC(id) = ic + if (ID_APS(id) != NULL) + ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1] + if (!override) + if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES) + next + + if (!trace) { + i = id_gid (id, start) + ID_LINE(id,1) = line[1] + ID_LINE(id,2) = line[2] + } + + if (IS_INDEFD(shift[1])) + fit_shift[1] = INDEFD + else { + if (!trace) + fit_shift[1] = fit_shift[1] - shift[1] + else + fit_shift[1] = -shift[1] + } + + ID_IC(id) = ic1 + call id_gdata (id) + iferr (call id_fitdata (id)) + ; + + call ri_loghdr (id, reference, logfd, nlogfd, loghdr) + loghdr = 0 + call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd) + + if (ID_NFEATURES(id) < nreid && trace) { + call ri_loghdr (id, reference, logfd, nlogfd, 3) + break + } + + if (ID_NFEATURES(id) > 0) { + call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO) + call id_saveid (id, line) + } + } + + ID_IC(id) = ic + i = id_gid (id, start) + fit_shift[1] = fit_shift[2] + for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1); + line[1]=line[1]+step[1]) { + ID_LINE(id,1) = line[1] + ID_AP(id,1) = line[1] + ID_IC(id) = ic + if (ID_APS(id) != NULL) + ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1] + if (!override) + if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES) + next + + if (!trace) { + i = id_gid (id, start) + ID_LINE(id,1) = line[1] + ID_LINE(id,2) = line[2] + } + + if (IS_INDEFD(shift[1])) + fit_shift[1] = INDEFD + else { + if (!trace) + fit_shift[1] = fit_shift[1] + shift[1] + else + fit_shift[1] = shift[1] + } + + ID_IC(id) = ic1 + call id_gdata (id) + iferr (call id_fitdata (id)) + ; + + call ri_loghdr (id, reference, logfd, nlogfd, loghdr) + loghdr = 0 + call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd) + + if (ID_NFEATURES(id) < nreid && trace) { + call ri_loghdr (id, reference, logfd, nlogfd, 3) + break + } + + if (ID_NFEATURES(id) > 0) { + call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO) + call id_saveid (id, line) + } + } + } + + + fit_shift[2] = 0. + for (line[2]=start[2]+step[2]; line[2]<=ID_MAXLINE(id,2); + line[2]=line[2]+step[2]) { + ID_LINE(id,2) = line[2] + ID_AP(id,2) = line[2] + ID_IC(id) = ic + + if (IS_INDEFD(shift[2])) + fit_shift[2] = INDEFD + else { + if (!trace) + fit_shift[2] = fit_shift[2] + shift[2] + else + fit_shift[2] = shift[2] + } + + fit_shift[1] = fit_shift[2] + for (line[1]=start[1]; line[1]>0; line[1]=line[1]-step[1]) { + ID_LINE(id,1) = line[1] + ID_AP(id,1) = line[1] + ID_IC(id) = ic + if (ID_APS(id) != NULL) + ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1] + if (!override) + if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES) + next + + if (!trace) { + i = id_gid (id, start) + ID_LINE(id,1) = line[1] + ID_LINE(id,2) = line[2] + } + + if (IS_INDEFD(shift[1])) + fit_shift[1] = INDEFD + else { + if (!trace) + fit_shift[1] = fit_shift[1] - shift[1] + else + fit_shift[1] = -shift[1] + } + + ID_IC(id) = ic1 + call id_gdata (id) + iferr (call id_fitdata (id)) + ; + + call ri_loghdr (id, reference, logfd, nlogfd, loghdr) + loghdr = 0 + call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd) + + if (ID_NFEATURES(id) < nreid && trace) { + call ri_loghdr (id, reference, logfd, nlogfd, 3) + break + } + + if (ID_NFEATURES(id) > 0) { + call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO) + call id_saveid (id, line) + } + } + + ID_IC(id) = ic + i = id_gid (id, start) + fit_shift[1] = fit_shift[2] + for (line[1]=start[1]+step[1]; line[1]<=ID_MAXLINE(id,1); + line[1]=line[1]+step[1]) { + ID_LINE(id,1) = line[1] + ID_AP(id,1) = line[1] + ID_IC(id) = ic + if (ID_APS(id) != NULL) + ID_AP(id,1) = Memi[ID_APS(id)+line[1]-1] + if (!override) + if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES) + next + + if (!trace) { + i = id_gid (id, start) + ID_LINE(id,1) = line[1] + ID_LINE(id,2) = line[2] + } + + if (IS_INDEFD(shift[1])) + fit_shift[1] = INDEFD + else { + if (!trace) + fit_shift[1] = fit_shift[1] + shift[1] + else + fit_shift[1] = shift[1] + } + + ID_IC(id) = ic1 + call id_gdata (id) + iferr (call id_fitdata (id)) + ; + + call ri_loghdr (id, reference, logfd, nlogfd, loghdr) + loghdr = 0 + call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd) + + if (ID_NFEATURES(id) < nreid && trace) { + call ri_loghdr (id, reference, logfd, nlogfd, 3) + break + } + + if (ID_NFEATURES(id) > 0) { + call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO) + call id_saveid (id, line) + } + } + } + + ID_IC(id) = ic + if (ic != ic1) + call ic_closed (ic1) + + call smw_close (MW(ID_SH(id))) + call imunmap (IM(ID_SH(id))) + call shdr_close (ID_SH(id)) +end + + +# RI_IMAGE -- Reidentify an image. + +procedure ri_image (id, reference, image, ans, logfd, nlogfd, pd) + +pointer id # ID pointer +char reference[ARB] # Reference image +char image[ARB] # Image to be reidentified +char ans[3] # Interactive? +int logfd[ARB] # Logfiles +int nlogfd # Number of logfiles +pointer pd # Plot file pointer + +bool newaps # Add new apertures not in reference? +bool override # Override previous identifications? +bool verbose # Verbose output? + +int i, j, ap, loghdr, id_getid(), id_dbcheck() +double shift, fit_shift, clgetd() +pointer ic, ic1 +bool clgetb() + +begin + # Open the image and return if there is an error. + call strcpy (image, Memc[ID_IMAGE(id)], SZ_FNAME) + iferr (call id_map (id)) { + call erract (EA_WARN) + return + } + call dtunmap (ID_DT(id)) + + newaps = clgetb ("newaps") + override = clgetb ("override") + verbose = clgetb ("verbose") + + ic = ID_IC(id) + if (ans[1] == 'N') + ic1 = ic + else + call ic_open (ic1) + + loghdr = 2 + shift = clgetd ("shift") + + # For MULTISPEC search the reference list of each aperture. If + # a reference of the same aperture is not found and the newaps + # flag is set use the initial reference and then add the + # reidentification to the reference list. + # For NDSPEC apply each reference to the image. + + if (SMW_FORMAT(MW(ID_SH(id))) == SMW_ES || + SMW_FORMAT(MW(ID_SH(id))) == SMW_MS) { + for (i=1; i<=ID_MAXLINE(id,1); i=i+1) { + ap = Memi[ID_APS(id)+i-1] + for (j=ID_NID(id); id_getid (id,j)!=EOF; j=j-1) + if (ap == ID_AP(id,1)) + break + + if (j == 0 && !newaps) { + if (verbose) { + call printf ( + "%s: Reference for aperture %d not found\n") + call pargstr (image) + call pargi (ap) + } + next + } + + ID_AP(id,1) = ap + ID_LINE(id,1) = i + + if (i == 1 && ic != ic1) + call ic_copy (ic, ic1) + + if (!override) + if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES) + next + + ID_IC(id) = ic1 + call id_gdata (id) + iferr (call id_fitdata (id)) + ; + + call ri_loghdr (id, reference, logfd, nlogfd, loghdr) + loghdr = 0 + + fit_shift = shift + call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd) + + if (ID_NFEATURES(id) > 0) { + call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO) + if (j == 0 && newaps) { + call id_sid (id, ID_NID(id)+1) + if (verbose) { + call printf ( + "%s: New reference for aperture %d\n") + call pargstr (image) + call pargi (ap) + } + } + } + ID_IC(id) = ic + } + + } else { + for (i=1; id_getid (id,i)!=EOF; i=i+1) { + if (i == 1 && ic != ic1) + call ic_copy (ic, ic1) + + if (!override) + if (id_dbcheck (id, Memc[ID_IMAGE(id)], ID_AP(id,1)) == YES) + next + + ID_IC(id) = ic1 + call id_gdata (id) + iferr (call id_fitdata (id)) + ; + + call ri_loghdr (id, reference, logfd, nlogfd, loghdr) + loghdr = 0 + + fit_shift = shift + call ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd) + + if (ID_NFEATURES(id) > 0) + call id_dbwrite (id, Memc[ID_IMAGE(id)], ID_AP(id,1), NO) + ID_IC(id) = ic + } + } + + ID_IC(id) = ic + if (ic != ic1) + call ic_closed (ic1) + call smw_close (MW(ID_SH(id))) + call imunmap (IM(ID_SH(id))) + call shdr_close (ID_SH(id)) +end + + +# RI_REIDENTIFY -- Reidentify features using a reference image database entry. + +procedure ri_reidentify (id, fit_shift, ans, logfd, nlogfd, pd) + +pointer id # ID pointer +double fit_shift # Shift in fit coords (input and output) +char ans[3] # Interactive? +int logfd[ARB] # Logfiles +int nlogfd # Number of logfiles +pointer pd # Plot file pointer + +int i, j, nfeatures1, nfeatures2, nfit, iden, mono, clgwrd() +double shift, pix_shift, z_shift, v, vrms +double id_fitpt(), fit_to_pix(), id_shift(), id_center(), id_rms(), id_zval() +pointer sp, pix, fit +bool clgetb() + +begin + nfeatures1 = ID_NFEATURES(id) + call smark (sp) + call salloc (pix, nfeatures1, TY_DOUBLE) + call salloc (fit, nfeatures1, TY_DOUBLE) + call amovd (PIX(id,1), Memd[pix], nfeatures1) + call amovd (FIT(id,1), Memd[fit], nfeatures1) + + # If no initial shift is given then the procedure id_shift + # computes a shift between the reference features and the + # features in the image. The purpose of the shift is to get the + # reference feature positions close enough to those of the image + # being identified that the centering algorithm will determine + # the exact positions of the features. An initial shift of zero + # is used if the two images are very nearly aligned as in the + # case of tracing features in a two dimensional image or for a + # set of images taken with the same observing setup. + + if (IS_INDEFD(fit_shift)) { + ID_FWIDTH(id) = FWIDTH(id,1) + ID_FTYPE(id) = FTYPE(id,1) + shift = id_shift (id) + } else + shift = fit_shift + + # For each reference feature a shift is added to bring the pixel + # position near that for the image being identified and then the + # centering algorithm is used. If the centering algorithm fails + # the feature is discarded. A mean shift is computed for the + # features which have been reidentified. + + do i = 1, ID_NFEATURES(id) { + PIX(id,i) = fit_to_pix (id, FIT(id,i) + shift) + PIX(id,i) = id_center (id, PIX(id,i), 1, FWIDTH(id,i), FTYPE(id,i), + NO) + if (!IS_INDEFD(PIX(id,i))) + FIT(id,i) = id_fitpt (id, PIX(id,i)) + } + for (i=1; i<ID_NFEATURES(id); i=i+1) { + if (IS_INDEFD(PIX(id,i))) + next + for (j=i+1; j<=ID_NFEATURES(id); j=j+1) { + if (IS_INDEFD(PIX(id,j))) + next + if (abs (PIX(id,i)-PIX(id,j)) < ID_MINSEP(id)) { + if (abs (FIT(id,i)-USER(id,i)) < abs (FIT(id,j)-USER(id,j))) + PIX(id,j) = INDEFD + else { + PIX(id,i) = INDEFD + break + } + } + } + } + + pix_shift = 0. + fit_shift = 0. + z_shift = 0. + j = 0 + do i = 1, ID_NFEATURES(id) { + if (IS_INDEFD(PIX(id,i))) + next + + pix_shift = pix_shift + PIX(id,i) - Memd[pix+i-1] + fit_shift = fit_shift + FIT(id,i) - Memd[fit+i-1] + if (Memd[fit+i-1] != 0.) + z_shift = z_shift + id_zval (id, FIT(id,i), Memd[fit+i-1]) + + j = j + 1 + PIX(id,j) = PIX(id,i) + FIT(id,j) = FIT(id,i) + USER(id,j) = USER(id,i) + WTS(id,j) = WTS(id,i) + FWIDTH(id,j) = FWIDTH(id,i) + FTYPE(id,j) = FTYPE(id,i) + } + ID_NFEATURES(id) = j + + nfeatures2 = j + pix_shift = pix_shift / max (1, ID_NFEATURES(id)) + fit_shift = fit_shift / max (1, ID_NFEATURES(id)) + z_shift = z_shift / max (1, ID_NFEATURES(id)) + + # If refitting the coordinate function is requested and there is + # more than one feature and there is a previously defined + # coordinate function then refit. Otherwise compute a coordinate + # shift. + + mono = YES + switch (ID_REFIT(id)) { + case 1: + if (ID_CV(id) != NULL) { + if (ID_NFEATURES(id)>1) + call id_dofit (id, NO) + else + call id_doshift (id, NO) + } + case 2: + if (ID_CV(id) != NULL) + call id_doshift (id, NO) + case 3: + call id_velocity (id, NO) + v = ID_REDSHIFT(id) * VLIGHT + vrms = ID_RMSRED(id) * VLIGHT + } + if (ID_NEWCV(id) == YES) { + iferr (call id_fitdata (id)) + mono = NO + call id_fitfeatures (id) + } + + if (clgetb ("addfeatures")) { + ID_FWIDTH(id) = FWIDTH(id,1) + ID_FTYPE(id) = FTYPE(id,1) + call id_linelist (id) + if (ID_NEWFEATURES(id) == YES) { + switch (ID_REFIT(id)) { + case 1: + if (ID_NFEATURES(id)>1) + call id_dofit (id, NO) + else + call id_doshift (id, NO) + case 2: + call id_doshift (id, NO) + case 3: + call id_velocity (id, NO) + v = ID_REDSHIFT(id) * VLIGHT + vrms = ID_RMSRED(id) * VLIGHT + } + if (ID_NEWCV(id) == YES) { + iferr (call id_fitdata (id)) + mono = NO + call id_fitfeatures (id) + } + } + } + + # Enter fitting interactively. + iden = NO + if ((ID_NFEATURES(id)>1) && (ID_CV(id)!=NULL || ID_REFIT(id) == 3)) { + if (ans[1] != 'N') { + if (ans[1] != 'Y') { + nfit = 0 + for (j=1; j<=ID_NFEATURES(id); j=j+1) + if (WTS(id,j) > 0.) + nfit = nfit + 1 + call printf ( + "%s%s%23t%3d/%-3d %3d/%-3d %9.3g %10.3g %7.3g %7.3g\n") + call pargstr (Memc[ID_IMAGE(id)]) + call pargstr (Memc[ID_SECTION(id)]) + call pargi (nfeatures2) + call pargi (nfeatures1) + call pargi (nfit) + call pargi (ID_NFEATURES(id)) + call pargd (pix_shift) + call pargd (fit_shift) + if (ID_REFIT(id) == 3) { + call pargd (v) + call pargd (vrms) + } else { + call pargd (z_shift) + call pargd (id_rms(id)) + } + call flush (STDOUT) + repeat { + ifnoerr (i = clgwrd ("answer", ans, SZ_FNAME, + "|no|yes|NO|YES|")) + break + } + call clpstr ("answer", ans) + } + switch (ans[1]) { + case 'y', 'Y': + mono = YES + i = ID_REFIT(id) + call reidentify (id) + ID_REFIT(id) = i + if (ID_REFIT(id) == 3) { + call id_velocity (id, NO) + v = ID_REDSHIFT(id) * VLIGHT + vrms = ID_RMSRED(id) * VLIGHT + } + iden = YES + } + if (ans[1] != 'Y') + call gdeactivate (ID_GP(id), 0) + } + } + + # Record log information if a log file descriptor is given. + for (i = 1; i <= nlogfd; i = i + 1) { + if (ans[1] == 'n' && logfd[i] == STDOUT) + next + nfit = 0 + for (j=1; j<=ID_NFEATURES(id); j=j+1) + if (WTS(id,j) > 0.) + nfit = nfit + 1 + call fprintf (logfd[i], + "%s%s%23t%3d/%-3d %3d/%-3d %9.3g %10.3g %7.3g %7.3g\n") + call pargstr (Memc[ID_IMAGE(id)]) + call pargstr (Memc[ID_SECTION(id)]) + call pargi (nfeatures2) + call pargi (nfeatures1) + call pargi (nfit) + call pargi (ID_NFEATURES(id)) + call pargd (pix_shift) + call pargd (fit_shift) + if (ID_REFIT(id) == 3) { + call pargd (v) + call pargd (vrms) + } else { + call pargd (z_shift) + call pargd (id_rms(id)) + } + if (ID_REFIT(id) == 3 && logfd[i] != STDOUT) + call id_log (id, "", logfd[i]) + if (mono == NO) + call fprintf (logfd[i], "Non-monotonic dispersion function") + call flush (logfd[i]) + if (logfd[i] == STDOUT) + iden = NO + } + # Print log if STDOUT is not used but if the IDENTIFY is done. + if (iden == YES) { + call printf ( + "%s%s%23t%3d/%-3d %3d/%-3d %9.3g %10.3g %7.3g %7.3g\n") + call pargstr (Memc[ID_IMAGE(id)]) + call pargstr (Memc[ID_SECTION(id)]) + call pargi (nfeatures2) + call pargi (nfeatures1) + call pargi (nfit) + call pargi (ID_NFEATURES(id)) + call pargd (pix_shift) + call pargd (fit_shift) + if (ID_REFIT(id) == 3) { + call pargd (v) + call pargd (vrms) + } else { + call pargd (z_shift) + call pargd (id_rms(id)) + } + if (mono == NO) + call printf ("Non-monotonic dispersion function") + call flush (STDOUT) + } + + # Make log plot. + call ri_plot (id, pd) + + call sfree (sp) +end + + +# RI_LOGHDR -- Print a log header in the log files. + +procedure ri_loghdr (id, reference, logfd, nlogfd, flag) + +pointer id # Identify structure +char reference[ARB] # Reference image +int logfd[ARB] # Log file descriptors +int nlogfd # Number of log files +int flag # Header type flag (1=banner, 2=Column labels, 3=Error) + +int i +pointer str + +begin + for (i = 1; i <= nlogfd; i = i + 1) { + switch (flag) { + case 1: # Print ID + call malloc (str, SZ_LINE, TY_CHAR) + call sysid (Memc[str], SZ_LINE) + switch (ID_REFIT(id)) { + case 1, 2: + call fprintf (logfd[i], "\nREIDENTIFY: %s\n") + call pargstr (Memc[str]) + case 3: + call fprintf (logfd[i], "\nRVREIDLINES: %s\n") + call pargstr (Memc[str]) + } + call mfree (str, TY_CHAR) + case 2: # Print labels + switch (ID_REFIT(id)) { + case 1, 2: + call fprintf (logfd[i], + " Reference image = %s, New image = %s, Refit = %b\n") + call pargstr (reference) + call pargstr (Memc[ID_IMAGE(id)]) + call pargi (ID_REFIT(id)) + call fprintf (logfd[i], + "%20s %7s %7s %9s %10s %7s %7s\n") + call pargstr ("Image Data") + call pargstr ("Found") + call pargstr ("Fit") + call pargstr ("Pix Shift") + call pargstr ("User Shift") + call pargstr ("Z Shift") + call pargstr ("RMS") + case 3: + call fprintf (logfd[i], + " Reference image = %s, New image = %s\n") + call pargstr (reference) + call pargstr (Memc[ID_IMAGE(id)]) + call fprintf (logfd[i], + "%20s %7s %7s %9s %10s %8s %7s\n") + call pargstr ("Image Data") + call pargstr ("Found") + call pargstr ("Fit") + call pargstr ("Pix Shift") + call pargstr ("User Shift") + call pargstr ("Velocity") + call pargstr ("RMS") + } + case 3: # Error + call fprintf (logfd[i], " ** Too many features lost **\n") + } + } +end + + +# RI_PLOT -- Plot residual graph of reidentified lines. + +procedure ri_plot (id, pd) + +pointer id # ID pointer +pointer pd # GIO pointer + +int i, j +pointer sp, str, x, y, gt, gt_init() +double id_zshiftd() + +begin + # Check if there is anything to plot. + if (pd == NULL || ID_NFEATURES(id) == 0) + return + + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (x, ID_NFEATURES(id), TY_REAL) + call salloc (y, ID_NFEATURES(id), TY_REAL) + + # Set plot points. + j = 0 + do i = 1, ID_NFEATURES(id) { + if (IS_INDEFD(USER(id,i))) + break + + Memr[x+j] = USER(id,i) + Memr[y+j] = id_zshiftd (id, FIT(id,i), 0) + j = j + 1 + } + + if (j == 0) { + call sfree (sp) + return + } + + # Make the plot. + call sprintf (Memc[str], SZ_LINE, "Reidentify: %s") + call pargstr (Memc[ID_IMAGE(id)]) + gt = gt_init () + call gt_sets (gt, GTTYPE, "mark") + call gt_sets (gt, GTXLABEL, "user coordinates") + call gt_sets (gt, GTYLABEL, "residuals (fit - user)") + call gt_sets (gt, GTTITLE, Memc[str]) + call gclear (pd) + call gascale (pd, Memr[x], j, 1) + call gascale (pd, Memr[y], j, 2) + call gt_swind (pd, gt) + call gt_labax (pd, gt) + call gt_plot (pd, gt, Memr[x], Memr[y], j) + call gt_free (gt) + + call sfree (sp) +end |