diff options
Diffstat (limited to 'noao/onedspec/ecidentify/t_ecreid.x')
| -rw-r--r-- | noao/onedspec/ecidentify/t_ecreid.x | 181 |
1 files changed, 181 insertions, 0 deletions
diff --git a/noao/onedspec/ecidentify/t_ecreid.x b/noao/onedspec/ecidentify/t_ecreid.x new file mode 100644 index 00000000..5e9769ff --- /dev/null +++ b/noao/onedspec/ecidentify/t_ecreid.x @@ -0,0 +1,181 @@ +include <error.h> +include <smw.h> +include "ecidentify.h" + +# T_ECREIDENTIFY -- Reidentify echelle features starting from reference. +# If no initial shift is given then the procedure ec_shift computes a +# shift between the reference features and the features in the image. +# The purpose of the shift is to get the feature positions from the +# reference image close enough to those of the image being identified +# that the centering algorithm will determine the exact positions of the +# features. The recentered features are then fit with either a shift +# of a full echelle function and written to database. + +procedure t_ecreidentify () + +int images # List of images +pointer ref # Reference image +double shift # Initial shift + +int i, j, fd, nfeatures1, nfeatures2 +double shift1, pix, fit, pix_shift, fit_shift, z_shift +pointer sp, log, ec + +int imtopenp(), ec_getim(), clpopnu(), clgfil(), open(), btoi() +double ec_fitpt(), ec_fittopix(), ec_shift(), ec_center(), ec_rms() +double clgetd() +bool clgetb() +real clgetr() +errchk ec_dbread(), ec_gdata(), ec_fitdata() + +begin + call smark (sp) + call salloc (ref, SZ_FNAME, TY_CHAR) + call salloc (log, SZ_FNAME, TY_CHAR) + + # Allocate the basic data structure. + call ec_init (ec) + + # Initialize fitting + call ecf_seti ("niterate", 0) + call ecf_setd ("low", 3.D0) + call ecf_setd ("high", 3.D0) + + # Get task parameters. + images = imtopenp ("images") + call clgstr ("reference", Memc[ref], SZ_FNAME) + shift = clgetd ("shift") + call clgstr ("database", Memc[EC_DATABASE(ec)], SZ_FNAME) + EC_CRADIUS(ec) = clgetr ("cradius") + EC_THRESHOLD(ec) = clgetr ("threshold") + EC_LOGFILES(ec) = clpopnu ("logfiles") + EC_REFIT(ec) = btoi (clgetb ("refit")) + + # Write logfile header. + while (clgfil (EC_LOGFILES(ec), Memc[log], SZ_FNAME) != EOF) { + iferr (fd = open (Memc[log], APPEND, TEXT_FILE)) { + call erract (EA_WARN) + next + } + call sysid (Memc[log], SZ_LINE) + call fprintf (fd, "\nECREIDENTIFY: %s\n") + call pargstr (Memc[log]) + call fprintf (fd, + " Reference image = %s, Refit = %b\n") + call pargstr (Memc[ref]) + call pargb (EC_REFIT(ec) == YES) + call fprintf (fd, "%20s %7s %7s %9s %10s %7s %7s\n") + call pargstr ("Image") + call pargstr ("Found") + call pargstr ("Fit") + call pargstr ("Pix Shift") + call pargstr ("User Shift") + call pargstr ("Z Shift") + call pargstr ("RMS") + call close (fd) + } + + # Reidentify features in each spectrum. + while (ec_getim (images, Memc[EC_IMAGE(ec)], SZ_FNAME) != EOF) { + call ec_gdata (ec) + call ec_dbread (ec, Memc[ref], NO) + call ec_fitdata (ec) + call ec_fitfeatures (ec) + + if (IS_INDEFD (shift)) { + EC_FWIDTH(ec) = FWIDTH(ec,1) + EC_FTYPE(ec) = abs (FTYPE(ec,1)) + EC_MINSEP(ec) = 1. + EC_MAXFEATURES(ec) = 20 + shift1 = ec_shift (ec) + } else + shift1 = shift + + # Recenter features. + pix_shift = 0. + fit_shift = 0. + z_shift = 0. + nfeatures1 = EC_NFEATURES(ec) + + j = 0. + do i = 1, EC_NFEATURES(ec) { + call ec_gline (ec, LINE(ec,i)) + pix = ec_fittopix (ec, FIT(ec,i) + shift1/ORDER(ec,i)) + pix = ec_center (ec, pix, FWIDTH(ec,i), FTYPE(ec,i)) + if (IS_INDEFD (pix)) + next + fit = ec_fitpt (ec, APN(ec,i), pix) + + pix_shift = pix_shift + pix - PIX(ec,i) + fit_shift = fit_shift + (fit - FIT(ec,i)) * ORDER(ec,i) + if (FIT(ec,i) != 0.) + z_shift = z_shift + (fit - FIT(ec,i)) / FIT(ec,i) + + j = j + 1 + APN(ec,j) = APN(ec,i) + LINE(ec,j) = LINE(ec,i) + ORDER(ec,j) = ORDER(ec,i) + PIX(ec,j) = pix + FIT(ec,j) = FIT(ec,i) + USER(ec,j) = USER(ec,i) + FWIDTH(ec,j) = FWIDTH(ec,i) + FTYPE(ec,j) = abs (FTYPE(ec,i)) + } + EC_NFEATURES(ec) = j + + # 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. + + if ((EC_REFIT(ec)==YES)&&(EC_NFEATURES(ec)>1)&&(EC_ECF(ec)!=NULL)) { + iferr (call ec_dofit (ec, NO, YES)) { + call erract (EA_WARN) + next + } + } else + call ec_doshift (ec, NO) + if (EC_NEWECF(ec) == YES) + call ec_fitfeatures (ec) + + nfeatures2 = 0 + do i = 1, EC_NFEATURES(ec) + if (FTYPE(ec,i) > 0) + nfeatures2 = nfeatures2 + 1 + + # Write a database entry for the reidentified image. + if (EC_NFEATURES(ec) > 0) + call ec_dbwrite (ec, Memc[EC_IMAGE(ec)], NO) + + # Record log information if a log file descriptor is given. + call clprew (EC_LOGFILES(ec)) + while (clgfil (EC_LOGFILES(ec), Memc[log], SZ_FNAME) != EOF) { + iferr (fd = open (Memc[log], APPEND, TEXT_FILE)) { + call erract (EA_WARN) + next + } + call fprintf (fd, + "%20s %3d/%-3d %3d/%-3d %9.3g %10.3g %7.3g %7.3g\n") + call pargstr (Memc[EC_IMAGE(ec)]) + call pargi (EC_NFEATURES(ec)) + call pargi (nfeatures1) + call pargi (nfeatures2) + call pargi (EC_NFEATURES(ec)) + call pargd (pix_shift / max (1, EC_NFEATURES(ec))) + call pargd (fit_shift / max (1, EC_NFEATURES(ec))) + call pargd (z_shift / max (1, EC_NFEATURES(ec))) + call pargd (ec_rms(ec)) + call close (fd) + } + + call smw_close (MW(EC_SH(ec))) + do i = 1, EC_NLINES(ec) + MW(SH(ec,i)) = NULL + } + + call dgsfree (EC_ECF(ec)) + call clpcls (EC_LOGFILES(ec)) + call ec_free (ec) + call imtclose (images) + call sfree (sp) +end |