aboutsummaryrefslogtreecommitdiff
path: root/noao/imred/crutil/src/crexamine.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/imred/crutil/src/crexamine.x')
-rw-r--r--noao/imred/crutil/src/crexamine.x626
1 files changed, 626 insertions, 0 deletions
diff --git a/noao/imred/crutil/src/crexamine.x b/noao/imred/crutil/src/crexamine.x
new file mode 100644
index 00000000..14e8d589
--- /dev/null
+++ b/noao/imred/crutil/src/crexamine.x
@@ -0,0 +1,626 @@
+include <error.h>
+include <syserr.h>
+include <imhdr.h>
+include <gset.h>
+include <mach.h>
+include <pkg/gtools.h>
+include "crlist.h"
+
+# CR_EXAMINE -- Examine cosmic ray candidates interactively.
+# CR_GRAPH -- Make a graph
+# CR_NEAREST -- Find the nearest cosmic ray to the cursor.
+# CR_DELETE -- Set replace flag for cosmic ray candidate nearest cursor.
+# CR_UNDELETE -- Set no replace flag for cosmic ray candidate nearest cursor.
+# CR_UPDATE -- Change replacement flags, thresholds, and graphs.
+# CR_PLOT -- Make log plot
+
+define HELP "crutil$src/cosmicrays.key"
+define PROMPT "cosmic ray options"
+
+# CR_EXAMINE -- Examine cosmic ray candidates interactively.
+
+procedure cr_examine (cr, gp, gt, im, fluxratio, first)
+
+pointer cr # Cosmic ray list
+pointer gp # GIO pointer
+pointer gt # GTOOLS pointer
+pointer im # Image pointer
+real fluxratio # Flux ratio threshold
+int first # Initial key
+
+char cmd[SZ_LINE]
+int i, newgraph, wcs, key, nc, nl, c1, c2, l1, l2, show
+real wx, wy, x1, y1, x2, y2
+pointer data
+
+int clgcur()
+pointer imgs2r()
+
+begin
+ # Set up the graphics.
+ call gt_sets (gt, GTPARAMS, IM_TITLE(im))
+
+ # Set image limits
+ nc = IM_LEN(im, 1)
+ nl = IM_LEN(im, 2)
+
+ # Enter cursor loop.
+ key = first
+ repeat {
+ switch (key) {
+ case '?': # Print help text.
+ call gpagefile (gp, HELP, PROMPT)
+ case ':': # Colon commands.
+ switch (cmd[1]) {
+ case '/':
+ call gt_colon (cmd, gp, gt, newgraph)
+ default:
+ call printf ("\007")
+ }
+ case 'a': # Toggle show all
+ if (show == 0)
+ show = 1
+ else
+ show = 0
+ newgraph = YES
+ case 'd': # Delete candidate
+ call cr_delete (gp, wx, wy, cr, i, show)
+ case 'e': # Delete candidates in region
+ x1 = wx; y1 = wy
+ call printf ("again:")
+ if (clgcur ("cursor", x2, y2, wcs, key, cmd, SZ_LINE) == EOF)
+ return
+ call cr_delete_reg (gp, x1, y1, x2, y2, cr, show)
+ case 'q': # Quit
+ break
+ case 'r': # Redraw the graph.
+ newgraph = YES
+ case 's': # Make surface plots
+ call cr_nearest (gp, wx, wy, cr, i, show)
+ c1 = max (1, int (Memr[CR_COL(cr)+i-1]) - 5)
+ c2 = min (nc, int (Memr[CR_COL(cr)+i-1]) + 5)
+ l1 = max (1, int (Memr[CR_LINE(cr)+i-1]) - 5)
+ l2 = min (nl, int (Memr[CR_LINE(cr)+i-1]) + 5)
+ data = imgs2r (im, c1, c2, l1, l2)
+ call gclear (gp)
+ call gsview (gp, 0.03, 0.48, 0.53, 0.98)
+ call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, -33., 25.)
+ call gsview (gp, 0.53, 0.98, 0.53, 0.98)
+ call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, -123., 25.)
+ call gsview (gp, 0.03, 0.48, 0.03, 0.48)
+ call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, 57., 25.)
+ call gsview (gp, 0.53, 0.98, 0.03, 0.48)
+ call cr_surface (gp, Memr[data], c2-c1+1, l2-l1+1, 147., 25.)
+ call fprintf (STDERR, "[Type any key to continue]")
+ i = clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE)
+ newgraph = YES
+ case 't': # Set threshold
+ call cr_update (gp, wy, cr, fluxratio, show)
+ call clputr ("fluxratio", fluxratio)
+ case 'u': # Undelete candidate
+ call cr_undelete (gp, wx, wy, cr, i, show)
+ case 'v': # Undelete candidates in region
+ x1 = wx; y1 = wy
+ call printf ("again:")
+ if (clgcur ("cursor", x2, y2, wcs, key, cmd, SZ_LINE) == EOF)
+ return
+ call cr_undelete_reg (gp, x1, y1, x2, y2, cr, show)
+ case 'w':# Window the graph.
+ call gt_window (gt, gp, "cursor", newgraph)
+ case ' ': # Print info
+ call cr_nearest (gp, wx, wy, cr, i, show)
+ call printf ("%d %d\n")
+ call pargr (Memr[CR_COL(cr)+i-1])
+ call pargr (Memr[CR_LINE(cr)+i-1])
+ case 'z': # NOP
+ newgraph = NO
+ default: # Ring bell for unrecognized commands.
+ call printf ("\007")
+ }
+
+ # Update the graph if needed.
+ if (newgraph == YES) {
+ call cr_graph (gp, gt, cr, fluxratio, show)
+ newgraph = NO
+ }
+ } until (clgcur ("cursor", wx, wy, wcs, key, cmd, SZ_LINE) == EOF)
+end
+
+
+# CR_GRAPH -- Make a graph
+
+procedure cr_graph (gp, gt, cr, fluxratio, show)
+
+pointer gp # GIO pointer
+pointer gt # GTOOLS pointers
+pointer cr # Cosmic ray list
+real fluxratio # Flux ratio threshold
+int show # Show (0=all, 1=train)
+
+int i, ncr
+real x1, x2, y1, y2
+pointer sp, x, y, w, flag, index
+
+begin
+ call smark (sp)
+
+ call cr_show (show, cr, x, y, w, flag, index, ncr)
+ if (ncr == 0) {
+ call sfree (sp)
+ return
+ }
+
+ call gclear (gp)
+ call gt_ascale (gp, gt, Memr[x+1], Memr[y+1], ncr)
+ call gt_swind (gp, gt)
+ call gt_labax (gp, gt)
+
+ do i = 1, ncr {
+ if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+ call gmark (gp, Memr[x+i], Memr[y+i], GM_PLUS, 2., 2.)
+ else
+ call gmark (gp, Memr[x+i], Memr[y+i], GM_CROSS, 2., 2.)
+ if (Memr[w+i] != 0.)
+ call gmark (gp, Memr[x+i], Memr[y+i], GM_BOX, 2., 2.)
+ }
+
+ call ggwind (gp, x1, x2, y1, y2)
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, x1, fluxratio, x2, fluxratio)
+
+ call sfree (sp)
+end
+
+
+# CR_NEAREST -- Find the nearest cosmic ray to the cursor.
+
+procedure cr_nearest (gp, wx, wy, cr, nearest, show)
+
+pointer gp # GIO pointer
+real wx, wy # Cursor position
+pointer cr # Cosmic ray list
+int nearest # Index of nearest point (returned)
+int show # Show (0=all, 1=train)
+
+int i, ncr
+real x0, y0, x1, y1, x2, y2, r2, r2min
+pointer sp, x, y, w, flag, index
+
+begin
+ call smark (sp)
+
+ call cr_show (show, cr, x, y, w, flag, index, ncr)
+ if (ncr == 0) {
+ call sfree (sp)
+ return
+ }
+
+ # Search for nearest point in NDC.
+ r2min = MAX_REAL
+ call gctran (gp, wx, wy, wx, wy, 1, 0)
+ do i = 1, ncr {
+ x1 = Memr[x+i]
+ y1 = Memr[y+i]
+ call gctran (gp, x1, y1, x0, y0, 1, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ x2 = x1
+ y2 = y1
+ nearest = i
+ }
+ }
+ if (index != NULL)
+ nearest = Memi[index+nearest]
+
+ # Move the cursor to the selected point.
+ call gscur (gp, x2, y2)
+
+ call sfree (sp)
+end
+
+
+# CR_DELETE -- Set replace flag for cosmic ray candidate nearest cursor.
+
+procedure cr_delete (gp, wx, wy, cr, nearest, show)
+
+pointer gp # GIO pointer
+real wx, wy # Cursor position
+pointer cr # Cosmic ray list
+int nearest # Index of nearest point (returned)
+int show # Show (0=all, 1=train)
+
+int i, ncr
+real x0, y0, x1, y1, x2, y2, r2, r2min
+pointer sp, x, y, w, flag, index
+
+begin
+ call smark (sp)
+
+ call cr_show (show, cr, x, y, w, flag, index, ncr)
+ if (ncr == 0) {
+ call sfree (sp)
+ return
+ }
+
+ # Search for nearest point in NDC.
+ nearest = 0
+ r2min = MAX_REAL
+ call gctran (gp, wx, wy, wx, wy, 1, 0)
+ do i = 1, ncr {
+ if ((Memi[flag+i] == YES) || (Memi[flag+i] == ALWAYSYES))
+ next
+ x1 = Memr[x+i]
+ y1 = Memr[y+i]
+ call gctran (gp, x1, y1, x0, y0, 1, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ x2 = x1
+ y2 = y1
+ nearest = i
+ }
+ }
+
+ # Move the cursor to the selected point and mark the deleted point.
+ if (nearest > 0) {
+ if (index != NULL)
+ nearest = Memi[index+nearest]
+ Memi[CR_FLAG(cr)+nearest-1] = ALWAYSYES
+ Memi[CR_WT(cr)+nearest-1] = -1
+ call gscur (gp, x2, y2)
+ call gseti (gp, G_PMLTYPE, 0)
+ y2 = Memr[CR_RATIO(cr)+nearest-1]
+ call gmark (gp, x2, y2, GM_PLUS, 2., 2.)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, x2, y2, GM_CROSS, 2., 2.)
+ }
+
+ call sfree (sp)
+end
+
+
+# CR_DELETE_REG -- Set replace flag for cosmic ray candidates in a region.
+
+procedure cr_delete_reg (gp, wx1, wy1, wx2, wy2, cr, show)
+
+pointer gp # GIO pointer
+real wx1, wy1, wx2, wy2 # Cursor positions
+pointer cr # Cosmic ray list
+int show # Show (0=all, 1=train)
+
+int i, j, ncr
+real x0, y0, x1, y1
+pointer sp, x, y, w, flag, index
+
+begin
+ call smark (sp)
+
+ call cr_show (show, cr, x, y, w, flag, index, ncr)
+ if (ncr == 0) {
+ call sfree (sp)
+ return
+ }
+
+ # Check order of region points.
+ if (wx1 > wx2) {
+ x0 = wx1
+ wx1 = wx2
+ wx2 = x0
+ }
+ if (wy1 > wy2) {
+ y0 = wy1
+ wy1 = wy2
+ wy2 = y0
+ }
+
+ # Check if point in region.
+ call gctran (gp, wx1, wy1, wx1, wy1, 1, 0)
+ call gctran (gp, wx2, wy2, wx2, wy2, 1, 0)
+ do i = 1, ncr {
+ if ((Memi[flag+i] == YES) || (Memi[flag+i] == ALWAYSYES))
+ next
+ x1 = Memr[x+i]
+ y1 = Memr[y+i]
+ call gctran (gp, x1, y1, x0, y0, 1, 0)
+
+ # Mark the deleted points.
+ if ((x0 > wx1) && (x0 < wx2) && (y0 > wy1) && (y0 < wy2)) {
+ if (index != NULL)
+ j = Memi[index+i]
+ else
+ j = i
+ Memi[CR_FLAG(cr)+j-1] = ALWAYSYES
+ Memi[CR_WT(cr)+j-1] = -1
+ call gscur (gp, x1, y1)
+ call gseti (gp, G_PMLTYPE, 0)
+ y1 = Memr[CR_RATIO(cr)+j-1]
+ call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+ }
+ }
+ call sfree (sp)
+end
+
+
+# CR_UNDELETE -- Set no replace flag for cosmic ray candidate nearest cursor.
+
+procedure cr_undelete (gp, wx, wy, cr, nearest, show)
+
+pointer gp # GIO pointer
+real wx, wy # Cursor position
+pointer cr # Cosmic ray list
+int nearest # Index of nearest point (returned)
+int show # Show (0=all, 1=train)
+
+int i, ncr
+real x0, y0, x1, y1, x2, y2, r2, r2min
+pointer sp, x, y, w, flag, index
+
+begin
+ call smark (sp)
+
+ call cr_show (show, cr, x, y, w, flag, index, ncr)
+ if (ncr == 0) {
+ call sfree (sp)
+ return
+ }
+
+ # Search for nearest point in NDC.
+ nearest = 0
+ r2min = MAX_REAL
+ call gctran (gp, wx, wy, wx, wy, 1, 0)
+ do i = 1, ncr {
+ if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+ next
+ x1 = Memr[x+i]
+ y1 = Memr[y+i]
+ call gctran (gp, x1, y1, x0, y0, 1, 0)
+ r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+ if (r2 < r2min) {
+ r2min = r2
+ x2 = x1
+ y2 = y1
+ nearest = i
+ }
+ }
+
+ # Move the cursor to the selected point and mark the delete point.
+ if (nearest > 0) {
+ if (index != NULL)
+ nearest = Memi[index+nearest]
+ Memi[CR_FLAG(cr)+nearest-1] = ALWAYSNO
+ Memi[CR_WT(cr)+nearest-1] = 1
+ call gscur (gp, x2, y2)
+
+ call gseti (gp, G_PMLTYPE, 0)
+ y2 = Memr[CR_RATIO(cr)+nearest-1]
+ call gmark (gp, x2, y2, GM_CROSS, 2., 2.)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, x2, y2, GM_PLUS, 2., 2.)
+ }
+
+ call sfree (sp)
+end
+
+
+# CR_UNDELETE_REG -- Set no replace flag for cosmic ray candidates in a region.
+
+procedure cr_undelete_reg (gp, wx1, wy1, wx2, wy2, cr, show)
+
+pointer gp # GIO pointer
+real wx1, wy1, wx2, wy2 # Cursor positions
+pointer cr # Cosmic ray list
+int show # Show (0=all, 1=train)
+
+int i, j, ncr
+real x0, y0, x1, y1
+pointer sp, x, y, w, flag, index
+
+begin
+ call smark (sp)
+
+ call cr_show (show, cr, x, y, w, flag, index, ncr)
+ if (ncr == 0) {
+ call sfree (sp)
+ return
+ }
+
+ # Check order of region points.
+ if (wx1 > wx2) {
+ x0 = wx1
+ wx1 = wx2
+ wx2 = x0
+ }
+ if (wy1 > wy2) {
+ y0 = wy1
+ wy1 = wy2
+ wy2 = y0
+ }
+
+ # Check if point in region.
+ call gctran (gp, wx1, wy1, wx1, wy1, 1, 0)
+ call gctran (gp, wx2, wy2, wx2, wy2, 1, 0)
+ do i = 1, ncr {
+ if ((Memi[flag+i] == NO) || (Memi[flag+i] == ALWAYSNO))
+ next
+ x1 = Memr[x+i]
+ y1 = Memr[y+i]
+ call gctran (gp, x1, y1, x0, y0, 1, 0)
+
+ # Mark the deleted points.
+ if ((x0 > wx1) && (x0 < wx2) && (y0 > wy1) && (y0 < wy2)) {
+ if (index != NULL)
+ j = Memi[index+i]
+ else
+ j = i
+ Memi[CR_FLAG(cr)+j-1] = ALWAYSNO
+ Memi[CR_WT(cr)+j-1] = 1
+ call gscur (gp, x1, y1)
+ call gseti (gp, G_PMLTYPE, 0)
+ y1 = Memr[CR_RATIO(cr)+j-1]
+ call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+ }
+ }
+ call sfree (sp)
+end
+
+
+# CR_UPDATE -- Change replacement flags, thresholds, and graphs.
+
+procedure cr_update (gp, wy, cr, fluxratio, show)
+
+pointer gp # GIO pointer
+real wy # Y cursor position
+pointer cr # Cosmic ray list
+real fluxratio # Flux ratio threshold
+int show # Show (0=all, 1=train)
+
+int i, ncr, flag
+real x1, x2, y1, y2
+pointer x, y, f
+
+begin
+ call gseti (gp, G_PLTYPE, 0)
+ call ggwind (gp, x1, x2, y1, y2)
+ call gline (gp, x1, fluxratio, x2, fluxratio)
+ fluxratio = wy
+ call gseti (gp, G_PLTYPE, 2)
+ call gline (gp, x1, fluxratio, x2, fluxratio)
+
+ if (show == 1)
+ return
+
+ ncr = CR_NCR(cr)
+ x = CR_FLUX(cr) - 1
+ y = CR_RATIO(cr) - 1
+ f = CR_FLAG(cr) - 1
+
+ do i = 1, ncr {
+ flag = Memi[f+i]
+ if ((flag == ALWAYSYES) || (flag == ALWAYSNO))
+ next
+ x1 = Memr[x+i]
+ y1 = Memr[y+i]
+ if (flag == NO) {
+ if (y1 < fluxratio) {
+ Memi[f+i] = YES
+ call gseti (gp, G_PMLTYPE, 0)
+ call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+ }
+ } else {
+ if (y1 >= fluxratio) {
+ Memi[f+i] = NO
+ call gseti (gp, G_PMLTYPE, 0)
+ call gmark (gp, x1, y1, GM_CROSS, 2., 2.)
+ call gseti (gp, G_PMLTYPE, 1)
+ call gmark (gp, x1, y1, GM_PLUS, 2., 2.)
+ }
+ }
+ }
+end
+
+
+# CR_PLOT -- Make log plot
+
+procedure cr_plot (cr, im, fluxratio)
+
+pointer cr # Cosmic ray list
+pointer im # Image pointer
+real fluxratio # Flux ratio threshold
+
+int fd, open(), errcode()
+pointer sp, fname, gp, gt, gopen(), gt_init()
+errchk gopen
+
+begin
+ call smark (sp)
+ call salloc (fname, SZ_FNAME, TY_CHAR)
+
+ # Open the plotfile.
+ call clgstr ("plotfile", Memc[fname], SZ_FNAME)
+ iferr (fd = open (Memc[fname], APPEND, BINARY_FILE)) {
+ if (errcode() != SYS_FNOFNAME)
+ call erract (EA_WARN)
+ return
+ }
+
+ # Set up the graphics.
+ gp = gopen ("stdplot", NEW_FILE, fd)
+ gt = gt_init()
+ call gt_sets (gt, GTTYPE, "mark")
+ call gt_sets (gt, GTXTRAN, "log")
+ call gt_setr (gt, GTXMIN, 10.)
+ call gt_setr (gt, GTYMIN, 0.)
+ call gt_sets (gt, GTTITLE, "Parameters of cosmic rays candidates")
+ call gt_sets (gt, GTPARAMS, IM_TITLE(im))
+ call gt_sets (gt, GTXLABEL, "Flux")
+ call gt_sets (gt, GTYLABEL, "Flux Ratio")
+
+ call cr_graph (gp, gt, cr, fluxratio, 'r')
+
+ call gt_free (gt)
+ call gclose (gp)
+ call close (fd)
+ call sfree (sp)
+end
+
+
+# CR_SHOW -- Select data to show.
+# This returns pointers to the data. Note the pointers are salloc from
+# the last smark which is done by the calling program.
+
+procedure cr_show (show, cr, x, y, w, flag, index, ncr)
+
+int show #I Data to show (0=all, 1=train)
+pointer cr #I CR data
+pointer x #O Fluxes
+pointer y #O Ratios
+pointer w #O Weights
+pointer flag #O Flags
+pointer index #O Index into CR data (if not null)
+int ncr #O Number of selected data points
+
+int i
+
+begin
+ switch (show) {
+ case 0:
+ ncr = CR_NCR(cr)
+ x = CR_FLUX(cr) - 1
+ y = CR_RATIO(cr) - 1
+ w = CR_WT(cr) - 1
+ flag = CR_FLAG(cr) - 1
+ index = NULL
+ case 1:
+ ncr = CR_NCR(cr)
+ call salloc (x, ncr, TY_REAL)
+ call salloc (y, ncr, TY_REAL)
+ call salloc (w, ncr, TY_REAL)
+ call salloc (flag, ncr, TY_INT)
+ call salloc (index, ncr, TY_INT)
+
+ ncr = 0
+ x = x - 1
+ y = y - 1
+ w = w - 1
+ flag = flag - 1
+ index = index - 1
+
+ do i = 1, CR_NCR(cr) {
+ if (Memr[CR_WT(cr)+i-1] == 0.)
+ next
+ ncr = ncr + 1
+ Memr[x+ncr] = Memr[CR_FLUX(cr)+i-1]
+ Memr[y+ncr] = Memr[CR_RATIO(cr)+i-1]
+ Memr[w+ncr] = Memr[CR_WT(cr)+i-1]
+ Memi[flag+ncr] = Memi[CR_FLAG(cr)+i-1]
+ Memi[index+ncr] = i
+ }
+ }
+end
generated by cgit (git 2.34.1) at 2025-09-20 09:08:29 -0400