diff options
Diffstat (limited to 'noao/imred/crutil/src')
24 files changed, 3842 insertions, 0 deletions
diff --git a/noao/imred/crutil/src/Revisions b/noao/imred/crutil/src/Revisions new file mode 100644 index 00000000..a5ef52f1 --- /dev/null +++ b/noao/imred/crutil/src/Revisions @@ -0,0 +1,151 @@ +.help revisions Nov99 crutil +.nf +t_cosmicrays.x + A pointer to a an array of pointers was used in one place as a real. This + is an error when integer and real arrays are not of the same size; i.e. + on 64-bit architectures. (8/2/12, Valdes) + +======= +V2.16.1 +======= + +t_craverage.x + The pointer for the input mask buffer when using the same mask for both + input and output was not being read to include the edge lines which are + used during the computation resulting in a segmentation violation. + (11/24/04, Valdes) + +======= +V2.12.2 +======= + +t_craverage.x + To grow an output mask create by the task requires closing the new mask + and reopening it READ_WRITE. (10/23/02, Valdes) + +======= +V2.12.1 +======= + +===== +V2.12 +===== + +t_craverage.x +t_crmedian.x + The mask name may be !<keyword>. (3/22/02, Valdes) + +xtmaskname.x + +t_craverage.x +t_crmedian.x +t_crgrow.x +mkpkg + Modified to allow FITS mask extensions. (3/22/02, Valdes) + +mkpkg + Added missing mkpkg depencies for several source files. (12/13/01, MJF) + +============================ +CRUTIL V1.5: August 22, 2001 +============================ + +crcombine.cl +crcombine.par +../doc/crcombine.hlp + Modified to use new version of IMCOMBINE. (8/22/01, Valdes) + +noao$imred/crutil/ + + Installed package into the NOAO.IMRED package. (8/22/01, Valdes) + +t_craverage.x + The nrej features was done wrong for nrej of 1 or 2. (5/3/01, Valdes) + +t_craverage.x + 1. The amov and aclr calls were using the wrong buffer start. + 2. Added missing imtclose calls. + 3. The growing needed to be moved outside the block buffering. + (10/4/00, Valdes as diagnosed by Davis) + +t_craverage.x +craverage.par +../doc/craverage.hlp + Added an nrej parameter to excluded additional pixels from the average. + This is needed to deal with cosmic rays which are bigger than one + pixel or very nearby additional cosmic rays. (9/13/00, Valdes) + +======================== +CRUTIL V1.4: Jan 6, 2000 +======================== + +t_crmedian.x + The calculation of the sigma would reference uninitialized data if + the image size was not a multiple of the sigma block size. + (1/6/00, Valdes) + +t_craverage.x + The indexing of the pixels to use for the sigma calculation was wrong. + (1/6/00, Valdes) + +t_craverage.x + +craverage.par + +x_crutil.e +mkpkg +../crutil.cl +../crutil.men +../crutil.hd + New task that finds cosmic rays against an average excluding the candidate + pixel. It also detects objects and prevents cosmic rays being detected + in them (i.e. the cores of stars). (11/30/99, Valdes) + +t_crgrow.x +crgrow.par + 1. Broke main grow loop into a subroutine that can be called by other + tasks. + 2. Added inval and outval parameters to allow selecting mask values. + (11/30/99, Valdes) + +======================== +CRUTIL V1.3: Oct 19, 1999 +======================== + +crnebula.cl + The rin and rout parameters were not being used and instead + the default values were hardwired. (10/19/99, Valdes) + +======================== +CRUTIL V1.2: Sep 4, 1998 +======================== + +t_crmedian.x + On images with more than about 500K pixels the median operation is + done in overlapping blocks of lines. The amount of overlap is + half of the line median size "lmed". The bug is that on output + the overlap regions end up being zero. The output i/o is now done + in non-overlapping blocks. (9/4/98, Valdes) + +===================== +CRUTIL V1.1: May 1998 +===================== + +crexamine.x +cosmicrays.key + +cosmicrays.hlp + The graphical deletion and undeletion of candidates now includes the + keys 'e' and 'v' to delete and undelete from a marked rectangular + region. Also the key file was moved to the source directory. + (Valdes/Shopbell 5/15/98) + +cosmicrays.par + Fixed type in "crmasks" parameter name. + (Valdes/Shopbell 5/15/98) + +=========== +CRUTIL V1.0 +=========== + +New package created April 24, 1998. + +======= +V2.11.1 +======= +.endhelp diff --git a/noao/imred/crutil/src/cosmicrays.key b/noao/imred/crutil/src/cosmicrays.key new file mode 100644 index 00000000..beac2835 --- /dev/null +++ b/noao/imred/crutil/src/cosmicrays.key @@ -0,0 +1,43 @@ + COSMIC RAY DETECTION AND REPLACEMENT + + +INTERACTIVE IMAGE CURSOR COMMANDS + +When using the image display cursor to define a list of training objects +the following keystrokes may be given. + +? Help +c Identify the object as a cosmic ray +s Identify the object as a star +g Switch to the graphics plot +q Quit and continue with the cleaning + +There are no colon commands. + + +INTERACTIVE GRAPHICS CURSOR COMMANDS + + The graph shows the ratio of the mean background subtracted flux within the +detection window (excluding the candidate cosmic ray and the second brightest +pixel) to the flux of the candidate cosmic ray pixel as a function of the +mean flux. Both coordinates have been multiplied by 100 so the ratio is in +precent. The peaks of extended objects have high flux ratios while true +cosmic rays have low ratios. The main purpose of this step is to set the +flux ratio threshold for discriminating stars and galaxies. The cursor +keys are: + +? Help +a Toggle between showing all objects and only the training objects +d Mark candidate for replacement (applys to '+' points) +e Mark candidates in a region for replacement (applys to '+' points) +q Quit. Returns to training or to replace the selected pixels +r Redraw the graph +s Make a surface plot for the candidate nearest the cursor +t Set the flux ratio threshold at the y cursor position +u Mark candidate to not be replaced (applys to 'x' points) +v Mark candidates in a region to not be replaced (applys to 'x' points) +w Adjust the graph window (see GTOOLS) +<space> Print the pixel coordinates for a candidate + +There are no colon commands except those for the windowing options (type +:\help or see GTOOLS). diff --git a/noao/imred/crutil/src/cosmicrays.par b/noao/imred/crutil/src/cosmicrays.par new file mode 100644 index 00000000..9016d9f9 --- /dev/null +++ b/noao/imred/crutil/src/cosmicrays.par @@ -0,0 +1,17 @@ +input,s,a,,,,List of images in which to detect cosmic rays +output,s,a,,,,List of cosmic ray replaced output images (optional) +crmasks,s,h,"",,,"List of bad pixel masks (optional) +" +threshold,r,h,25.,,,Detection threshold above mean +fluxratio,r,h,2.,,,Flux ratio threshold (in percent) +npasses,i,h,5,1,,Number of detection passes +window,s,h,"5","5|7",,"Size of detection window +" +interactive,b,h,yes,,,Examine parameters interactively? +train,b,h,no,,,Use training objects? +objects,*imcur,h,"",,,Cursor list of training objects +savefile,f,h,"",,,File to save train objects +plotfile,f,h,"",,,Plot file +graphics,s,h,"stdgraph",,,Interactive graphics output device +cursor,*gcur,h,"",,,Graphics cursor input +answer,s,q,,"no|yes|NO|YES",,Review parameters for a particular image? diff --git a/noao/imred/crutil/src/craverage.par b/noao/imred/crutil/src/craverage.par new file mode 100644 index 00000000..b67c4585 --- /dev/null +++ b/noao/imred/crutil/src/craverage.par @@ -0,0 +1,23 @@ +input,f,a,,,,List of input images +output,f,a,,,,List of output images +crmask,f,h,,,,List of output cosmic ray and object masks +average,f,h,,,,List of output block average filtered images +sigma,f,h,,,,"List of output sigma images +" +navg,i,h,5,3,,Block average box size +nrej,i,h,0,0,,Number of high pixels to reject from the average +nbkg,i,h,5,1,,Background annulus width +nsig,i,h,25,10,,Box size for sigma calculation +var0,r,h,0.,0.,,Variance coefficient for DN^0 term +var1,r,h,0.,0.,,Variance coefficient for DN^1 term +var2,r,h,0.,0.,,"Variance coefficient for DN^2 term +" +crval,i,h,1,0,,Mask value for cosmic rays +lcrsig,r,h,10.,,,Low cosmic ray sigma outside object +hcrsig,r,h,5.,,,High cosmic ray sigma outside object +crgrow,r,h,0.,0.,,"Cosmic ray grow radius +" +objval,i,h,0,0,,Mask value for objects +lobjsig,r,h,10.,,,Low object detection sigma +hobjsig,r,h,5.,,,High object detection sigma +objgrow,r,h,0.,0.,,Object grow radius diff --git a/noao/imred/crutil/src/crcombine.cl b/noao/imred/crutil/src/crcombine.cl new file mode 100644 index 00000000..7c6bdc77 --- /dev/null +++ b/noao/imred/crutil/src/crcombine.cl @@ -0,0 +1,17 @@ +# CRCOMBINE -- Reject cosmic rays by combining multiple exposures. + +procedure crcombine (input, output) + +begin + imcombine (input, output, logfile=logfile, combine=combine, + reject=reject, scale=scale, zero=zero, statsec=statsec, + lsigma=lsigma, hsigma=hsigma, rdnoise=rdnoise, gain=gain, + grow=grow, headers=headers, bpmasks=bpmasks, rejmasks=rejmasks, + nrejmasks=nrejmasks, expmasks=expmasks, sigmas=sigmas, + project=project, outtype=outtype, outlimits=outlimits, + offsets=offsets, masktype=masktype, maskvalue=maskvalue, + blank=blank, weight=weight, expname=expname, + lthreshold=lthreshold, hthreshold=hthreshold, nlow=nlow, + nhigh=nhigh, nkeep=nkeep, mclip=mclip, snoise=snoise, + sigscale=sigscale, pclip=pclip) +end diff --git a/noao/imred/crutil/src/crcombine.par b/noao/imred/crutil/src/crcombine.par new file mode 100644 index 00000000..95ae95ae --- /dev/null +++ b/noao/imred/crutil/src/crcombine.par @@ -0,0 +1,45 @@ +# CRCOMBINE -- Image combine parameters + +input,s,a,,,,List of images to combine +output,s,a,,,,List of output images +logfile,s,h,"STDOUT",,,"Log file + +# Cosmic ray rejection parameters" +combine,s,h,"average","average|median|sum",,Type of combine operation +reject,s,h,"crreject","none|minmax|ccdclip|crreject|sigclip|avsigclip|pclip",,Type of rejection +scale,s,h,"mode",,,Image scaling +zero,s,h,"none",,,Image zero point offset +statsec,s,h,"",,,Image section for computing statistics +lsigma,r,h,10.,0.,,Lower sigma clipping factor +hsigma,r,h,3.,0.,,Upper sigma clipping factor +rdnoise,s,h,"0.",,,CCD readout noise (electrons) +gain,s,h,"1.",,,CCD gain (electrons/DN) +grow,r,h,0.,0.,,"Radius (pixels) for neighbor rejection + +# Additional output" +headers,s,h,"",,,List of header files (optional) +bpmasks,s,h,"",,,List of bad pixel masks (optional) +rejmasks,s,h,"",,,List of rejection masks (optional) +nrejmasks,s,h,"",,,List of number rejected masks (optional) +expmasks,s,h,"",,,List of exposure masks (optional) +sigmas,s,h,"",,,"List of sigma images (optional) + +# Additional parameters" +project,b,h,no,,,Project highest dimension of input images? +outtype,s,h,"real","short|ushort|integer|long|real|double",,Output image pixel datatype +outlimits,s,h,"",,,Output limits (x1 x2 y1 y2 ...) +offsets,f,h,"none",,,Input image offsets +masktype,s,h,"none","none|goodvalue|badvalue|goodbits|badbits",,Mask type +maskvalue,r,h,0,,,Mask value +blank,r,h,0.,,,Value if there are no pixels +weight,s,h,"none",,,Image weights +expname,s,h,"",,,Image header exposure time keyword +lthreshold,r,h,INDEF,,,Lower threshold +hthreshold,r,h,INDEF,,,Upper threshold +nlow,i,h,1,0,,minmax: Number of low pixels to reject +nhigh,i,h,1,0,,minmax: Number of high pixels to reject +nkeep,i,h,1,,,Minimum to keep (pos) or maximum to reject (neg) +mclip,b,h,yes,,,Use median in sigma clipping algorithms? +snoise,s,h,"0.",,,Sensitivity noise (fraction) +sigscale,r,h,0.1,0.,,Tolerance for sigma clipping scaling corrections +pclip,r,h,-0.5,,,pclip: Percentile clipping parameter diff --git a/noao/imred/crutil/src/credit.cl b/noao/imred/crutil/src/credit.cl new file mode 100644 index 00000000..8b844f33 --- /dev/null +++ b/noao/imred/crutil/src/credit.cl @@ -0,0 +1,13 @@ +# CREDIT -- Edit cosmic rays with an image display. + +procedure credit (input, output) + +begin + imedit (input, output, cursor=cursor, logfile=logfile, + display=display, autodisplay=autodisplay, + autosurface=autosurface, aperture=aperture, radius=radius, + search=search, buffer=buffer, width=width, xorder=xorder, + yorder=yorder, value=value, sigma=sigma, angh=angh, angv=angv, + command=command, graphics=graphics, default=default, + fixpix=fixpix) +end diff --git a/noao/imred/crutil/src/credit.par b/noao/imred/crutil/src/credit.par new file mode 100644 index 00000000..a23e0d35 --- /dev/null +++ b/noao/imred/crutil/src/credit.par @@ -0,0 +1,22 @@ +input,s,a,,,,Images to be edited +output,s,a,,,,Output images +cursor,*imcur,h,"",,,Cursor input +logfile,s,h,"",,,Logfile for record of cursor commands +display,b,h,yes,,,Display images? +autodisplay,b,h,yes,,,Automatic image display? +autosurface,b,h,no,,,Automatic surface plots? +aperture,s,h,"circular","|circular|square|",,Aperture type +radius,r,h,2.,,,Substitution radius +search,r,h,2.,,,Search radius +buffer,r,h,1.,0.,,Background buffer width +width,r,h,2.,1.,,Background width +xorder,i,h,2,0,,Background x order +yorder,i,h,2,0,,Background y order +value,r,h,0.,,,Constant value substitution +sigma,r,h,INDEF,,,Added noise sigma +angh,r,h, -33.,,,Horizontal viewing angle (degrees) +angv,r,h,25.,,,Vertical viewing angle (degrees) +command,s,h,"display $image 1 erase=$erase fill=yes order=0 >& dev$null",,,Display command +graphics,s,h,"stdgraph",,,Graphics device +default,s,h,"b",,,Default option for x-y input +fixpix,b,h,no,,,Fixpix style input? 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 diff --git a/noao/imred/crutil/src/crfind.x b/noao/imred/crutil/src/crfind.x new file mode 100644 index 00000000..58850940 --- /dev/null +++ b/noao/imred/crutil/src/crfind.x @@ -0,0 +1,305 @@ +include <math/gsurfit.h> + +# CR_FIND -- Find cosmic ray candidates. +# This procedure is an interface to special procedures specific to a given +# window size. + +procedure cr_find (cr, threshold, data, nc, nl, col, line, + sf1, sf2, x, y, z, w) + +pointer cr # Cosmic ray list +real threshold # Detection threshold +pointer data[ARB] # Data lines +int nc # Number of columns +int nl # Number of lines +int col # First column +int line # Center line +pointer sf1, sf2 # Surface fitting +real x[ARB], y[ARB], z[ARB], w[ARB] # Surface arrays + +pointer a, b, c, d, e, f, g + +begin + switch (nl) { + case 5: + a = data[1] + b = data[2] + c = data[3] + d = data[4] + e = data[5] + call cr_find5 (cr, threshold, col, line, Memr[a], Memr[b], + Memr[c], Memr[d], Memr[e], nc, sf1, sf2, x, y, z, w) + case 7: + a = data[1] + b = data[2] + c = data[3] + d = data[4] + e = data[5] + f = data[6] + g = data[7] + call cr_find7 (cr, threshold, col, line, Memr[a], Memr[b], + Memr[c], Memr[d], Memr[e], Memr[f], Memr[g], nc, + sf1, sf2, x, y, z, w) + } +end + + +# CR_FIND7 -- Find cosmic rays candidates in 7x7 window. +# This routine finds cosmic rays candidates with the following algorithm. +# 1. If the pixel is not a local maximum relative to it's 48 neighbors +# go on to the next pixel. +# 2. Identify the next strongest pixel in the 7x7 region. +# This suspect pixel is excluded in the following. +# 2. Compute the flux of the 7x7 region excluding the cosmic ray +# candidate and the suspect pixel. +# 3. The candidate must exceed the average flux per pixel by a specified +# threshold. If not go on to the next pixel. +# 4. Fit a plane to the border pixels (excluding the suspect pixel). +# 5. Subtract the background defined by the plane. +# 6. Determine a replacement value as the average of the four adjacent +# pixels (excluding the suspect pixels). +# 7. Add the pixel to the cosmic ray candidate list. + +procedure cr_find7 (cr, threshold, col, line, a, b, c, d, e, f, g, n, + sf1, sf2, x, y, z, w) + +pointer cr # Cosmic ray list +real threshold # Detection threshold +int col # First column +int line # Line +real a[ARB], b[ARB], c[ARB], d[ARB] # Image lines +real e[ARB], f[ARB], g[ARB] # Image lines +int n # Number of columns +pointer sf1, sf2 # Surface fitting +real x[49], y[49], z[49], w[49] # Surface arrays + +real bkgd[49] +int i1, i2, i3, i4, i5, i6, i7, j, j1, j2 +real p, flux, replace, asumr() +pointer sf + +begin + for (i4=4; i4<=n-3; i4=i4+1) { + # Must be local maxima. + p = d[i4] + if (p<a[i4]||p<b[i4]||p<c[i4]||p<e[i4]||p<f[i4]||p<g[i4]) + next + i1 = i4 - 3 + if (p<a[i1]||p<b[i1]||p<c[i1]||p<d[i1]||p<e[i1]||p<f[i1]||p<g[i1]) + next + i2 = i4 - 2 + if (p<a[i2]||p<b[i2]||p<c[i2]||p<d[i2]||p<e[i2]||p<f[i2]||p<g[i2]) + next + i3 = i4 - 1 + if (p<a[i3]||p<b[i3]||p<c[i3]||p<d[i3]||p<e[i3]||p<f[i3]||p<g[i3]) + next + i5 = i4 + 1 + if (p<a[i5]||p<b[i5]||p<c[i5]||p<d[i5]||p<e[i5]||p<f[i5]||p<g[i5]) + next + i6 = i4 + 2 + if (p<a[i6]||p<b[i6]||p<c[i6]||p<d[i6]||p<e[i6]||p<f[i6]||p<g[i6]) + next + i7 = i4 + 3 + if (p<a[i7]||p<b[i7]||p<c[i7]||p<d[i7]||p<e[i7]||p<f[i7]||p<g[i7]) + next + + # Convert to a single array in surface fitting order. + call amovr (a[i1], z[1], 7) + z[8] = b[i7]; z[9] = c[i7]; z[10] = d[i7]; z[11] = e[i7] + z[12] = f[i7]; z[13] = g[i7]; z[14] = g[i6]; z[15] = g[i5] + z[16] = f[i4]; z[17] = g[i3]; z[18] = g[i2]; z[19] = g[i1] + z[20] = f[i1]; z[21] = e[i1]; z[22] = d[i1]; z[23] = c[i1] + z[24] = b[i1] + call amovr (b[i2], z[25], 5) + call amovr (c[i2], z[30], 5) + call amovr (d[i2], z[35], 5) + call amovr (e[i2], z[40], 5) + call amovr (f[i2], z[45], 5) + + # Find the highest point excluding the center. + j1 = 37; j2 = 1 + do j = 2, 49 { + if (j == j1) + next + if (z[j] > z[j2]) + j2 = j + } + + # Compute the flux excluding the extreme points. + flux = (asumr (z, 49) - z[j1] - z[j2]) / 47 + + # Pixel must be exceed specified threshold. + if (p < flux + threshold) + next + + # Fit and subtract the background. + if (j2 < 25) { + w[j2] = 0 + sf = sf2 + call gsfit (sf, x, y, z, w, 24, WTS_USER, j) + w[j2] = 1 + } else { + sf = sf1 + call gsrefit (sf, x, y, z, w, j) + } + + call gsvector (sf, x, y, bkgd, 49) + call asubr (z, bkgd, z, 49) + p = z[j1] + + # Compute the flux excluding the extreme points. + flux = (asumr (z, 49) - z[j1] - z[j2]) / 47 + + # Determine replacement value from four nearest neighbors again + # excluding the most deviant pixels. + replace = 0 + j = 0 + if (j2 != 32) { + replace = replace + c[i4] + j = j + 1 + } + if (j2 != 36) { + replace = replace + d[i3] + j = j + 1 + } + if (j2 != 38) { + replace = replace + d[i5] + j = j + 1 + } + if (j2 != 42) { + replace = replace + e[i4] + j = j + 1 + } + replace = replace / j + + # Add pixel to cosmic ray list. + flux = 100. * flux + call cr_add (cr, col+i4-1, line, flux, flux/p, 0., replace, 0) + i4 = i7 + } +end + + +# CR_FIND5 -- Find cosmic rays candidates in 5x5 window. +# This routine finds cosmic rays candidates with the following algorithm. +# 1. If the pixel is not a local maximum relative to it's 24 neighbors +# go on to the next pixel. +# 2. Identify the next strongest pixel in the 5x5 region. +# This suspect pixel is excluded in the following. +# 2. Compute the flux of the 5x5 region excluding the cosmic ray +# candidate and the suspect pixel. +# 3. The candidate must exceed the average flux per pixel by a specified +# threshold. If not go on to the next pixel. +# 4. Fit a plane to the border pixels (excluding the suspect pixel). +# 5. Subtract the background defined by the plane. +# 6. Determine a replacement value as the average of the four adjacent +# pixels (excluding the suspect pixels). +# 7. Add the pixel to the cosmic ray candidate list. + +procedure cr_find5 (cr, threshold, col, line, a, b, c, d, e, n, + sf1, sf2, x, y, z, w) + +pointer cr # Cosmic ray list +real threshold # Detection threshold +int col # First column +int line # Line +real a[ARB], b[ARB], c[ARB], d[ARB], e[ARB] # Image lines +int n # Number of columns +pointer sf1, sf2 # Surface fitting +real x[25], y[25], z[25], w[25] # Surface arrays + +real bkgd[25] +int i1, i2, i3, i4, i5, j, j1, j2 +real p, flux, replace, asumr() +pointer sf + +begin + for (i3=3; i3<=n-2; i3=i3+1) { + # Must be local maxima. + p = c[i3] + if (p<a[i3]||p<b[i3]||p<d[i3]||p<e[i3]) + next + i1 = i3 - 2 + if (p<a[i1]||p<b[i1]||p<c[i1]||p<d[i1]||p<e[i1]) + next + i2 = i3 - 1 + if (p<a[i2]||p<b[i2]||p<c[i2]||p<d[i2]||p<e[i2]) + next + i4 = i3 + 1 + if (p<a[i4]||p<b[i4]||p<c[i4]||p<d[i4]||p<e[i4]) + next + i5 = i3 + 2 + if (p<a[i5]||p<b[i5]||p<c[i5]||p<d[i5]||p<e[i5]) + next + + # Convert to a single array in surface fitting order. + call amovr (a[i1], z[1], 5) + z[6] = b[i5]; z[7] = c[i5]; z[8] = d[i5]; z[9] = e[i5] + z[10] = e[i4]; z[11] = e[i3]; z[12] = e[i2]; z[13] = e[i1] + z[14] = d[i1]; z[15] = c[i1]; z[16] = b[i1] + call amovr (b[i2], z[17], 3) + call amovr (c[i2], z[20], 3) + call amovr (d[i2], z[23], 3) + + # Find the highest point excluding the center. + j1 = 21; j2 = 1 + do j = 2, 25 { + if (j == j1) + next + if (z[j] > z[j2]) + j2 = j + } + + # Compute the flux excluding the extreme points. + flux = (asumr (z, 25) - z[j1] - z[j2]) / 23 + + # Pixel must be exceed specified threshold. + if (p < flux + threshold) + next + + # Fit and subtract the background. + if (j2 < 17) { + w[j2] = 0 + sf = sf2 + call gsfit (sf, x, y, z, w, 16, WTS_USER, j) + w[j2] = 1 + } else { + sf = sf1 + call gsrefit (sf, x, y, z, w, j) + } + + call gsvector (sf, x, y, bkgd, 25) + call asubr (z, bkgd, z, 25) + p = z[j1] + + # Compute the flux excluding the extreme points. + flux = (asumr (z, 25) - z[j1] - z[j2]) / 23 + + # Determine replacement value from four nearest neighbors again + # excluding the most deviant pixels. + replace = 0 + j = 0 + if (j2 != 18) { + replace = replace + b[i3] + j = j + 1 + } + if (j2 != 20) { + replace = replace + c[i2] + j = j + 1 + } + if (j2 != 22) { + replace = replace + c[i4] + j = j + 1 + } + if (j2 != 24) { + replace = replace + d[i3] + j = j + 1 + } + replace = replace / j + + # Add pixel to cosmic ray list. + flux = 100. * flux + call cr_add (cr, col+i3-1, line, flux, flux/p, 0., replace, 0) + i3 = i5 + } +end diff --git a/noao/imred/crutil/src/crfix.cl b/noao/imred/crutil/src/crfix.cl new file mode 100644 index 00000000..68947c7a --- /dev/null +++ b/noao/imred/crutil/src/crfix.cl @@ -0,0 +1,20 @@ +# CRFIX -- Replace cosmic rays in an image using a cosmic ray mask. + +procedure crfix (input, output, crmask) + +file input {prompt="Input image"} +file output {prompt="Output image"} +file crmask {prompt="Cosmic ray mask"} + +begin + file in, out, crm + + in = input + out = output + crm = crmask + + if (in != out) + imcopy (in, out, verbose=no) + fixpix (out, crm, linterp="INDEF", cinterp="INDEF", verbose=no, + pixels=no) +end diff --git a/noao/imred/crutil/src/crgrow.par b/noao/imred/crutil/src/crgrow.par new file mode 100644 index 00000000..f57eb3cc --- /dev/null +++ b/noao/imred/crutil/src/crgrow.par @@ -0,0 +1,7 @@ +# CRGROW + +input,s,a,,,,Input cosmic ray masks +output,s,a,,,,Output cosmic ray masks +radius,r,h,1.,,,Grow radius +inval,i,h,INDEF,,,Input mask value to grow (INDEF for any) +outval,i,h,INDEF,,,Output grown mask value (INDEF for any) diff --git a/noao/imred/crutil/src/crlist.h b/noao/imred/crutil/src/crlist.h new file mode 100644 index 00000000..1ed498a7 --- /dev/null +++ b/noao/imred/crutil/src/crlist.h @@ -0,0 +1,17 @@ +define CR_ALLOC 100 # Allocation block size +define CR_LENSTRUCT 9 # Length of structure + +define CR_NCR Memi[$1] # Number of cosmic rays +define CR_NALLOC Memi[$1+1] # Length of cosmic ray list +define CR_COL Memi[$1+2] # Pointer to columns +define CR_LINE Memi[$1+3] # Pointer to lines +define CR_FLUX Memi[$1+4] # Pointer to fluxes +define CR_RATIO Memi[$1+5] # Pointer to flux ratios +define CR_WT Memi[$1+6] # Pointer to training weights +define CR_REPLACE Memi[$1+7] # Pointer to replacement values +define CR_FLAG Memi[$1+8] # Pointer to rejection flag + +define ALWAYSNO 3 +define ALWAYSYES 4 + +define CR_RMAX 3. # Maximum radius for matching diff --git a/noao/imred/crutil/src/crlist.x b/noao/imred/crutil/src/crlist.x new file mode 100644 index 00000000..bb49fb03 --- /dev/null +++ b/noao/imred/crutil/src/crlist.x @@ -0,0 +1,417 @@ +include <error.h> +include <syserr.h> +include <gset.h> +include <pmset.h> +include "crlist.h" + +define HELP "noao$lib/scr/cosmicrays.key" +define PROMPT "cosmic ray options" + +# CR_OPEN -- Open cosmic ray list +# CR_CLOSE -- Close cosmic ray list +# CR_ADD -- Add a cosmic ray candidate to cosmic ray list. +# CR_TRAIN -- Set flux ratio threshold from a training set. +# CR_FINDTHRESH -- Find flux ratio. +# CR_WEIGHT -- Compute the training weight at a particular flux ratio. +# CR_FLAGS -- Set cosmic ray reject flags. +# CR_BADPIX -- Store cosmic rays in bad pixel list. +# CR_REPLACE -- Replace cosmic rays in image with replacement values. + +# CR_OPEN -- Open cosmic ray list + +procedure cr_open (cr) + +pointer cr # Cosmic ray list pointer +errchk malloc + +begin + call malloc (cr, CR_LENSTRUCT, TY_STRUCT) + call malloc (CR_COL(cr), CR_ALLOC, TY_REAL) + call malloc (CR_LINE(cr), CR_ALLOC, TY_REAL) + call malloc (CR_FLUX(cr), CR_ALLOC, TY_REAL) + call malloc (CR_RATIO(cr), CR_ALLOC, TY_REAL) + call malloc (CR_WT(cr), CR_ALLOC, TY_REAL) + call malloc (CR_REPLACE(cr), CR_ALLOC, TY_REAL) + call malloc (CR_FLAG(cr), CR_ALLOC, TY_INT) + CR_NCR(cr) = 0 + CR_NALLOC(cr) = CR_ALLOC +end + + +# CR_CLOSE -- Close cosmic ray list + +procedure cr_close (cr) + +pointer cr # Cosmic ray list pointer + +begin + call mfree (CR_COL(cr), TY_REAL) + call mfree (CR_LINE(cr), TY_REAL) + call mfree (CR_FLUX(cr), TY_REAL) + call mfree (CR_RATIO(cr), TY_REAL) + call mfree (CR_WT(cr), TY_REAL) + call mfree (CR_REPLACE(cr), TY_REAL) + call mfree (CR_FLAG(cr), TY_INT) + call mfree (cr, TY_STRUCT) +end + +# CR_ADD -- Add a cosmic ray candidate to cosmic ray list. + +procedure cr_add (cr, col, line, flux, ratio, wt, replace, flag) + +pointer cr # Cosmic ray list pointer +int col # Cofluxn +int line # Line +real flux # Luminosity +real ratio # Ratio +real wt # Weight +real replace # Sky value +int flag # Flag value + +int ncr +errchk realloc + +begin + if (CR_NCR(cr) == CR_NALLOC(cr)) { + CR_NALLOC(cr) = CR_NALLOC(cr) + CR_ALLOC + call realloc (CR_COL(cr), CR_NALLOC(cr), TY_REAL) + call realloc (CR_LINE(cr), CR_NALLOC(cr), TY_REAL) + call realloc (CR_FLUX(cr), CR_NALLOC(cr), TY_REAL) + call realloc (CR_RATIO(cr), CR_NALLOC(cr), TY_REAL) + call realloc (CR_WT(cr), CR_NALLOC(cr), TY_REAL) + call realloc (CR_REPLACE(cr), CR_NALLOC(cr), TY_REAL) + call realloc (CR_FLAG(cr), CR_NALLOC(cr), TY_INT) + } + + ncr = CR_NCR(cr) + CR_NCR(cr) = ncr + 1 + Memr[CR_COL(cr)+ncr] = col + Memr[CR_LINE(cr)+ncr] = line + Memr[CR_FLUX(cr)+ncr] = flux + Memr[CR_RATIO(cr)+ncr] = ratio + Memr[CR_WT(cr)+ncr] = wt + Memr[CR_REPLACE(cr)+ncr] = replace + Memi[CR_FLAG(cr)+ncr] = flag +end + + +# CR_TRAIN -- Set flux ratio threshold from a training set. + +procedure cr_train (cr, gp, gt, im, fluxratio, fname) + +pointer cr #I Cosmic ray list +pointer gp #I GIO pointer +pointer gt #I GTOOLS pointer +pointer im #I IMIO pointer +real fluxratio #O Flux ratio threshold +char fname[ARB] #I Save file name + +char cmd[10] +bool gflag +real x, y, y1, y2, w, r, rmin +int i, j, n, f, ncr, wcs, key, fd, clgcur(), open(), errcode() +pointer col, line, ratio, flux, wt, flag + +begin + # Open save file + iferr (fd = open (fname, APPEND, TEXT_FILE)) { + if (errcode() != SYS_FNOFNAME) + call erract (EA_WARN) + fd = 0 + } + + ncr = CR_NCR(cr) + col = CR_COL(cr) - 1 + line = CR_LINE(cr) - 1 + flux = CR_FLUX(cr) - 1 + ratio = CR_RATIO(cr) - 1 + wt = CR_WT(cr) - 1 + flag = CR_FLAG(cr) - 1 + + gflag = false + n = 0 + while (clgcur ("objects", x, y, wcs, key, cmd, 10) != EOF) { + switch (key) { + case '?': + call gpagefile (gp, HELP, PROMPT) + next + case 'q': + break + case 's': + w = 1 + f = ALWAYSNO + case 'c': + w = -1 + f = ALWAYSYES + case 'g': + if (gflag) + call cr_examine (cr, gp, gt, im, fluxratio, 'z') + else { + if (n > 1) + call cr_findthresh (cr, fluxratio) + call cr_flags (cr, fluxratio) + call cr_examine (cr, gp, gt, im, fluxratio, 'r') + gflag = true + } + next + default: + next + } + + y1 = y - CR_RMAX + y2 = y + CR_RMAX + for (i=10; i<ncr && y1>Memr[line+i]; i=i+10) + ; + j = i - 9 + rmin = (Memr[col+j] - x) ** 2 + (Memr[line+j] - y) ** 2 + for (i=j+1; i<ncr && y2>Memr[line+i]; i=i+1) { + r = (Memr[col+i] - x) ** 2 + (Memr[line+i] - y) ** 2 + if (r < rmin) { + rmin = r + j = i + } + } + if (sqrt (rmin) > CR_RMAX) + next + + Memr[wt+j] = w + Memi[flag+j] = f + n = n + 1 + + if (gflag) { + if (n > 1) { + call cr_findthresh (cr, r) + call cr_update (gp, r, cr, fluxratio, 0) + } + call gmark (gp, Memr[flux+j], Memr[ratio+j], GM_BOX, 2., 2.) + } + if (fd > 0) { + call fprintf (fd, "%g %g %d %c\n") + call pargr (x) + call pargr (y) + call pargi (wcs) + call pargi (key) + } + } + + if (fd > 0) + call close (fd) +end + + +# CR_FINDTHRESH -- Find flux ratio. + +procedure cr_findthresh (cr, fluxratio) + +pointer cr #I Cosmic ray list +real fluxratio #O Flux ratio threshold + +real w, r, rmin, cr_weight() +int i, ncr +pointer ratio, wt + +begin + ncr = CR_NCR(cr) + ratio = CR_RATIO(cr) - 1 + wt = CR_WT(cr) - 1 + + fluxratio = Memr[ratio+1] + rmin = cr_weight (fluxratio, Memr[ratio+1], Memr[wt+1], ncr) + do i = 2, ncr { + if (Memr[wt+i] == 0.) + next + r = Memr[ratio+i] + w = cr_weight (r, Memr[ratio+1], Memr[wt+1], ncr) + if (w <= rmin) { + if (w == rmin) + fluxratio = min (fluxratio, r) + else { + rmin = w + fluxratio = r + } + } + } +end + + +# CR_WEIGHT -- Compute the training weight at a particular flux ratio. + +real procedure cr_weight (fluxratio, ratio, wts, ncr) + +real fluxratio #I Flux ratio +real ratio[ARB] #I Ratio Values +real wts[ARB] #I Weights +int ncr #I Number of ratio values +real wt #O Sum of weights + +int i + +begin + wt = 0. + do i = 1, ncr { + if (ratio[i] > fluxratio) { + if (wts[i] < 0.) + wt = wt - wts[i] + } else { + if (wts[i] > 0.) + wt = wt + wts[i] + } + } + return (wt) +end + + +# CR_FLAGS -- Set cosmic ray reject flags. + +procedure cr_flags (cr, fluxratio) + +pointer cr # Cosmic ray candidate list +real fluxratio # Rejection limits + +int i, ncr +pointer ratio, flag + +begin + ncr = CR_NCR(cr) + ratio = CR_RATIO(cr) - 1 + flag = CR_FLAG(cr) - 1 + + do i = 1, ncr { + if ((Memi[flag+i] == ALWAYSYES) || (Memi[flag+i] == ALWAYSNO)) + next + if (Memr[ratio+i] > fluxratio) + Memi[flag+i] = NO + else + Memi[flag+i] = YES + } +end + + +# CR_BADPIX -- Store cosmic rays in bad pixel list. +# This is currently a temporary measure until a real bad pixel list is +# implemented. + +procedure cr_badpix (cr, fname) + +pointer cr # Cosmic ray list +char fname[ARB] # Bad pixel file name + +int i, ncr, c, l, f, fd, open(), errcode() +pointer col, line, ratio, flux, flag +errchk open + +begin + # Open bad pixel file + iferr (fd = open (fname, APPEND, TEXT_FILE)) { + if (errcode() != SYS_FNOFNAME) + call erract (EA_WARN) + return + } + + ncr = CR_NCR(cr) + col = CR_COL(cr) - 1 + line = CR_LINE(cr) - 1 + flux = CR_FLUX(cr) - 1 + ratio = CR_RATIO(cr) - 1 + flag = CR_FLAG(cr) - 1 + + do i = 1, ncr { + f = Memi[flag+i] + if ((f == NO) || (f == ALWAYSNO)) + next + + c = Memr[col+i] + l = Memr[line+i] + call fprintf (fd, "%d %d\n") + call pargi (c) + call pargi (l) + } + call close (fd) +end + + +# CR_CRMASK -- Store cosmic rays in cosmic ray mask. + +procedure cr_crmask (cr, fname, ref) + +pointer cr # Cosmic ray list +char fname[ARB] # Cosmic ray mask +pointer ref # Refrence image + +int i, axlen[7], depth, ncr, f, rl[3,2] +long v[2] +pointer sp, title, pm, col, line, flag, pm_newmask() +errchk pm_loadf, pm_newmask, pm_savef + +begin + call smark (sp) + call salloc (title, SZ_LINE, TY_CHAR) + + pm = pm_newmask (ref, 1) + iferr (call pm_loadf (pm, fname, Memc[title], SZ_LINE)) + ; + call pm_gsize (pm, i, axlen, depth) + + ncr = CR_NCR(cr) + col = CR_COL(cr) - 1 + line = CR_LINE(cr) - 1 + flag = CR_FLAG(cr) - 1 + + RL_LEN(rl) = 2 + RL_AXLEN(rl) = axlen[1] + RL_N(rl,2) = 1 + RL_V(rl,2) = 1 + RL_X(rl,2) = 1 + + do i = 1, ncr { + f = Memi[flag+i] + if ((f == NO) || (f == ALWAYSNO)) + next + + v[1] = Memr[col+i] + v[2] = Memr[line+i] + call pmplri (pm, v, rl, depth, 1, PIX_SRC) + } + + call strcpy ("Cosmic ray mask from COSMICRAYS", Memc[title], SZ_LINE) + call pm_savef (pm, fname, Memc[title], PM_UPDATE) + call pm_close (pm) + call sfree (sp) +end + + +# CR_REPLACE -- Replace cosmic rays in image with replacement values. + +procedure cr_replace (cr, offset, im, nreplaced) + +pointer cr # Cosmic ray list +int offset # Offset in list +pointer im # IMIO pointer of output image +int nreplaced # Number replaced (for log) + +int i, ncr, c, l, f +real r +pointer col, line, replace, flag, imps2r() + +begin + ncr = CR_NCR(cr) + if (ncr <= offset) + return + + col = CR_COL(cr) - 1 + line = CR_LINE(cr) - 1 + replace = CR_REPLACE(cr) - 1 + flag = CR_FLAG(cr) - 1 + + do i = offset+1, ncr { + f = Memi[flag+i] + if ((f == NO) || (f == ALWAYSNO)) + next + + c = Memr[col+i] + l = Memr[line+i] + r = Memr[replace+i] + Memr[imps2r (im, c, c, l, l)] = r + nreplaced = nreplaced + 1 + } +end diff --git a/noao/imred/crutil/src/crmedian.par b/noao/imred/crutil/src/crmedian.par new file mode 100644 index 00000000..2baccb05 --- /dev/null +++ b/noao/imred/crutil/src/crmedian.par @@ -0,0 +1,15 @@ +input,f,a,,,,Input image +output,f,a,,,,Output image +crmask,f,h,,,,Output cosmic ray mask +median,f,h,,,,Output median image +sigma,f,h,,,,Output sigma image +residual,f,h,,,,Output residual image +var0,r,h,0.,0.,,Variance coefficient for DN^0 term +var1,r,h,0.,0.,,Variance coefficient for DN^1 term +var2,r,h,0.,0.,,Variance coefficient for DN^2 term +lsigma,r,h,10.,,,Low clipping sigma factor +hsigma,r,h,3.,,,High clipping sigma factor +ncmed,i,h,5,1,,Column box size for median level calculation +nlmed,i,h,5,1,,Line box size for median level calculation +ncsig,i,h,25,10,,Column box size for sigma calculation +nlsig,i,h,25,10,,Line box size for sigma calculation diff --git a/noao/imred/crutil/src/crnebula.cl b/noao/imred/crutil/src/crnebula.cl new file mode 100644 index 00000000..c5f214f9 --- /dev/null +++ b/noao/imred/crutil/src/crnebula.cl @@ -0,0 +1,116 @@ +# CRNEBULA -- Cosmic ray cleaning for images with fine nebular structure. + +procedure crnebula (input, output) + +file input {prompt="Input image"} +file output {prompt="Output image"} +file crmask {prompt="Cosmic ray mask"} +file residual {prompt="Residual median image"} +file rmedresid {prompt="Residual ring median image"} +real var0 = 1. {prompt="Variance coefficient for DN^0 term", min=0.} +real var1 = 0. {prompt="Variance coefficient for DN^1 term", min=0.} +real var2 = 0. {prompt="Variance coefficient for DN^2 term", min=0.} +real sigmed = 3 {prompt="Sigma clip factor for residual median"} +real sigrmed = 3 {prompt="Sigma clip factor for residual ring median"} + +int mbox = 5 {prompt="Median box size"} +real rin = 1.5 {prompt="Inner radius for ring median"} +real rout = 6 {prompt="Outer radius for ring median"} +bool verbose = no {prompt="Verbose"} + +begin + file in, out, med, sig, resid, rmed + struct expr + + # Query once for query parameters. + in = input + out = output + + # Check on output images. + if (out != "" && imaccess (out)) + error (1, "Output image already exists ("//out//")") + if (crmask != "" && imaccess (crmask)) + error (1, "Output mask already exists ("//crmask//")") + if (residual != "" && imaccess (residual)) + error (1, "Output residual image already exists ("//residual//")") + if (rmedresid != "" && imaccess (rmedresid)) + error (1, + "Output ring median difference already exists ("//rmedresid//")") + + # Create median results. + med = mktemp ("cr") + sig = mktemp ("cr") + resid = residual + if (resid == "") + resid = mktemp ("cr") + if (verbose) + printf ("Creating CRMEDIAN results\n") + crmedian (in, "", crmask="", median=med, sigma=sig, residual=resid, + var0=var0, var1=var1, var2=var2, lsigma=100., hsigma=sigmed, + ncmed=mbox, nlmed=mbox, ncsig=25, nlsig=25) + + # Create ring median filtered image. + rmed = mktemp ("cr") + rmedian (in, rmed, rin, rout, ratio=1., theta=0., zloreject=INDEF, + zhireject=INDEF, boundary="wrap", constant=0., verbose=verbose) + + # Create output images. + if (rmedresid != "") { + printf ("(a-b)/c\n") | scan (expr) + imexpr (expr, rmedresid, med, rmed, sig, dims="auto", + intype="auto", outtype="real", refim="auto", bwidth=0, + btype="nearest", bpixval=0., rangecheck=yes, verbose=no, + exprdb="none") + imdelete (rmed, verify-) + imdelete (sig, verify-) + if (out != "") { + if (verbose) + printf ("Create output image %s\n", out) + printf ("((a<%.3g)||(abs(b)>%.3g)) ? c : d\n", + sigmed, sigrmed) | scan (expr) + imexpr (expr, out, resid, rmedresid, in, med, dims="auto", + intype="auto", outtype="auto", refim="auto", bwidth=0, + btype="nearest", bpixval=0., rangecheck=yes, verbose=no, + exprdb="none") + } + if (crmask != "") { + if (verbose) + printf ("Create cosmic ray mask %s\n", crmask) + printf ("((a<%.3g)||(abs(b)>%.3g)) ? 0 : 1\n", + sigmed, sigrmed) | scan (expr) + set imtype = "pl" + imexpr (expr, crmask, resid, rmedresid, dims="auto", + intype="auto", outtype="short", refim="auto", bwidth=0, + btype="nearest", bpixval=0., rangecheck=yes, verbose=no, + exprdb="none") + } + imdelete (med, verify-) + } else { + if (out != "") { + if (verbose) + printf ("Create output image %s\n", out) + printf ("((a<%.3g)||(abs((b-c)/d)>%.3g)) ? e : b\n", + sigmed, sigrmed) | scan (expr) + imexpr (expr, out, resid, med, rmed, sig, in, dims="auto", + intype="auto", outtype="auto", refim="auto", bwidth=0, + btype="nearest", bpixval=0., rangecheck=yes, verbose=no, + exprdb="none") + } + if (crmask != "") { + if (verbose) + printf ("Create cosmic ray mask %s\n", crmask) + printf ("((a<%.3g)||(abs((b-c)/d)>%.3g)) ? 0 : 1\n", + sigmed, sigrmed) | scan (expr) + set imtype = "pl" + imexpr (expr, crmask, resid, med, rmed, sig, dims="auto", + intype="auto", outtype="short", refim="auto", bwidth=0, + btype="nearest", bpixval=0., rangecheck=yes, verbose=no, + exprdb="none") + } + imdelete (med, verify-) + imdelete (sig, verify-) + imdelete (rmed, verify-) + } + if (residual == "") + imdelete (resid, verify-) +end diff --git a/noao/imred/crutil/src/crsurface.x b/noao/imred/crutil/src/crsurface.x new file mode 100644 index 00000000..32645ff4 --- /dev/null +++ b/noao/imred/crutil/src/crsurface.x @@ -0,0 +1,46 @@ +define DUMMY 6 + +# CR_SURFACE -- Draw a perspective view of a surface. The altitude +# and azimuth of the viewing angle are variable. + +procedure cr_surface(gp, data, ncols, nlines, angh, angv) + +pointer gp # GIO pointer +real data[ncols,nlines] # Surface data to be plotted +int ncols, nlines # Dimensions of surface +real angh, angv # Orientation of surface (degrees) + +int wkid +pointer sp, work + +int first +real vpx1, vpx2, vpy1, vpy2 +common /frstfg/ first +common /noaovp/ vpx1, vpx2, vpy1, vpy2 + +begin + call smark (sp) + call salloc (work, 2 * (2 * ncols * nlines + ncols + nlines), TY_REAL) + + # Initialize surface common blocks + first = 1 + call srfabd() + + # Define viewport. + call ggview (gp, vpx1, vpx2, vpy1, vpy2) + + # Link GKS to GIO + wkid = 1 + call gopks (STDERR) + call gopwk (wkid, DUMMY, gp) + call gacwk (wkid) + + call ezsrfc (data, ncols, nlines, angh, angv, Memr[work]) + + call gdawk (wkid) + # We don't want to close the GIO pointer. + #call gclwk (wkid) + call gclks () + + call sfree (sp) +end diff --git a/noao/imred/crutil/src/mkpkg b/noao/imred/crutil/src/mkpkg new file mode 100644 index 00000000..1279b656 --- /dev/null +++ b/noao/imred/crutil/src/mkpkg @@ -0,0 +1,38 @@ +# COSMIC RAY CLEANING PACKAGE + +$call relink +$exit + +update: + $call relink + $call install + ; + +relink: + $update libpkg.a + $call crutil + ; + +install: + $move xx_crutil.e noaobin$x_crutil.e + ; + +crutil: + $omake x_crutil.x + $link x_crutil.o libpkg.a -lxtools -lcurfit -lgsurfit -lncar -lgks\ + -o xx_crutil.e + ; + +libpkg.a: + crexamine.x crlist.h <error.h> <gset.h> <mach.h> <pkg/gtools.h>\ + <imhdr.h> + crfind.x <math/gsurfit.h> + crlist.x crlist.h <error.h> <gset.h> <pmset.h> + crsurface.x + t_cosmicrays.x crlist.h <error.h> <math/gsurfit.h> <imhdr.h> <gset.h>\ + <pkg/gtools.h> <imset.h> + t_craverage.x <imhdr.h> <error.h> <mach.h> + t_crgrow.x <error.h> <imhdr.h> + t_crmedian.x <imhdr.h> <mach.h> + xtmaskname.x + ; diff --git a/noao/imred/crutil/src/t_cosmicrays.x b/noao/imred/crutil/src/t_cosmicrays.x new file mode 100644 index 00000000..fa68a39d --- /dev/null +++ b/noao/imred/crutil/src/t_cosmicrays.x @@ -0,0 +1,329 @@ +include <error.h> +include <imhdr.h> +include <imset.h> +include <math/gsurfit.h> +include <gset.h> +include <pkg/gtools.h> +include "crlist.h" + +# T_COSMICRAYS -- Detect and remove cosmic rays in images. +# A list of images is examined for cosmic rays which are then replaced +# by values from neighboring pixels. The output image may be the same +# as the input image. This is the top level procedure which manages +# the input and output image data. The actual algorithm for detecting +# cosmic rays is in CR_FIND. + +procedure t_cosmicrays () + +int list1 # List of input images to be cleaned +int list2 # List of output images +int list3 # List of output bad pixel files +real threshold # Detection threshold +real fluxratio # Luminosity boundary for stars +int npasses # Number of cleaning passes +int szwin # Size of detection window +bool train # Use training objects? +pointer savefile # Save file for training objects +bool interactive # Examine cosmic ray parameters? +char ans # Answer to interactive query + +int nc, nl, c, c1, c2, l, l1, l2, szhwin, szwin2 +int i, j, k, m, ncr, ncrlast, nreplaced, flag +pointer sp, input, output, badpix, str, gp, gt, im, in, out +pointer x, y, z, w, sf1, sf2, cr, data, ptr + +bool clgetb(), streq(), strne() +char clgetc() +int imtopenp(), imtlen(), imtgetim(), clpopnu(), clgfil(), clgeti() +real clgetr() +pointer immap(), impl2r(), imgs2r(), gopen(), gt_init() +errchk immap, impl2r, imgs2r +errchk cr_find, cr_examine, cr_replace, cr_plot, cr_crmask + +begin + call smark (sp) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (badpix, SZ_FNAME, TY_CHAR) + call salloc (savefile, SZ_FNAME, TY_CHAR) + call salloc (str, SZ_LINE, TY_CHAR) + + # Get the task parameters. Check that the number of output images + # is either zero, in which case the cosmic rays will be removed + # in place, or equal to the number of input images. + + list1 = imtopenp ("input") + list2 = imtopenp ("output") + i = imtlen (list1) + j = imtlen (list2) + if (j > 0 && j != i) + call error (0, "Input and output image lists do not match") + + list3 = clpopnu ("crmasks") + threshold = clgetr ("threshold") + fluxratio = clgetr ("fluxratio") + npasses = clgeti ("npasses") + szwin = clgeti ("window") + train = clgetb ("train") + call clgstr ("savefile", Memc[savefile], SZ_FNAME) + interactive = clgetb ("interactive") + call clpstr ("answer", "yes") + ans = 'y' + + # Set up the graphics. + call clgstr ("graphics", Memc[str], SZ_LINE) + if (interactive) { + gp = gopen (Memc[str], NEW_FILE+AW_DEFER, STDGRAPH) + 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, GTXLABEL, "Flux") + call gt_sets (gt, GTYLABEL, "Flux Ratio") + } + + # Set up surface fitting. The background points are placed together + # at the beginning of the arrays. There are two surface pointers, + # one for using the fast refit if there are no points excluded and + # one for doing a full fit with points excluded. + + szhwin = szwin / 2 + szwin2 = szwin * szwin + call salloc (data, szwin, TY_INT) + call salloc (x, szwin2, TY_REAL) + call salloc (y, szwin2, TY_REAL) + call salloc (z, szwin2, TY_REAL) + call salloc (w, szwin2, TY_REAL) + + k = 0 + do i = 1, szwin { + Memr[x+k] = i + Memr[y+k] = 1 + k = k + 1 + } + do i = 2, szwin { + Memr[x+k] = szwin + Memr[y+k] = i + k = k + 1 + } + do i = szwin-1, 1, -1 { + Memr[x+k] = i + Memr[y+k] = szwin + k = k + 1 + } + do i = szwin-1, 2, -1 { + Memr[x+k] = 1 + Memr[y+k] = i + k = k + 1 + } + do i = 2, szwin-1 { + do j = 2, szwin-1 { + Memr[x+k] = j + Memr[y+k] = i + k = k + 1 + } + } + call aclrr (Memr[z], szwin2) + call amovkr (1., Memr[w], 4*szwin-4) + call gsinit (sf1, GS_POLYNOMIAL, 2, 2, NO, 1., real(szwin), + 1., real(szwin)) + call gsinit (sf2, GS_POLYNOMIAL, 2, 2, NO, 1., real(szwin), + 1., real(szwin)) + call gsfit (sf1, Memr[x], Memr[y], Memr[z], Memr[w], 4*szwin-4, + WTS_USER, j) + + # Process each input image. Either work in place or create a + # new output image. If an error mapping the images occurs + # issue a warning and go on to the next input image. + + while (imtgetim (list1, Memc[input], SZ_FNAME) != EOF) { + if (imtgetim (list2, Memc[output], SZ_FNAME) == EOF) + call strcpy (Memc[input], Memc[output], SZ_FNAME) + if (clgfil (list3, Memc[badpix], SZ_FNAME) == EOF) + Memc[badpix] = EOS + + iferr { + in = NULL + out = NULL + cr = NULL + + # Map the input image. If the output image is + # the same as the input image work in place. + # Initialize IMIO to use a scrolling buffer of lines. + + if (streq (Memc[input], Memc[output])) { + im = immap (Memc[input], READ_WRITE, 0) + } else + im = immap (Memc[input], READ_ONLY, 0) + in = im + + nc = IM_LEN(in,1) + nl = IM_LEN(in,2) + if ((nl < szwin) || (nc < szwin)) + call error (0, "Image size is too small") + call imseti (in, IM_NBUFS, szwin) + call imseti (in, IM_TYBNDRY, BT_NEAREST) + call imseti (in, IM_NBNDRYPIX, szhwin) + + # Open the output image if needed. + if (strne (Memc[input], Memc[output])) + im = immap (Memc[output], NEW_COPY, in) + out = im + + # Open a cosmic ray list structure. + call cr_open (cr) + ncrlast = 0 + nreplaced = 0 + + # Now proceed through the image line by line, scrolling + # the line buffers at each step. If creating a new image + # also write out each line as it is read. A procedure is + # called to find the cosmic ray candidates in the line + # and add them to the list maintained by CRLIST. + # Note that cosmic rays are not replaced at this point + # in order to allow the user to modify the criteria for + # a cosmic ray and review the results. + + c1 = 1-szhwin + c2 = nc+szhwin + do i = 1, szwin-1 + Memi[data+i] = + imgs2r (in, c1, c2, i-szhwin, i-szhwin) + + do l = 1, nl { + do i = 1, szwin-1 + Memi[data+i-1] = Memi[data+i] + Memi[data+szwin-1] = + imgs2r (in, c1, c2, l+szhwin, l+szhwin) + if (out != in) + call amovr (Memr[Memi[data+szhwin]+szhwin], + Memr[impl2r(out,l)], nc) + + call cr_find (cr, threshold, Memi[data], + c2-c1+1, szwin, c1, l, + sf1, sf2, Memr[x], Memr[y], Memr[z], Memr[w]) + } + if (interactive && train) { + call cr_train (cr, gp, gt, in, fluxratio, Memc[savefile]) + train = false + } + call cr_flags (cr, fluxratio) + + # If desired examine the cosmic ray list interactively. + if (interactive && ans != 'N') { + if (ans != 'Y') { + call eprintf ("%s - ") + call pargstr (Memc[input]) + call flush (STDERR) + ans = clgetc ("answer") + } + if ((ans == 'Y') || (ans == 'y')) + call cr_examine (cr, gp, gt, in, fluxratio, 'r') + } + + # Now replace the selected cosmic rays in the output image. + + call imflush (out) + call imseti (out, IM_ADVICE, RANDOM) + call cr_replace (cr, ncrlast, out, nreplaced) + + # Do additional passes through the data. We work in place + # in the output image. Note that we only have to look in + # the vicinity of replaced cosmic rays for secondary + # events since we've already looked at every pixel once. + # Instead of scrolling through the image we will extract + # subrasters around each replaced cosmic ray. However, + # we use pointers into the subraster to maintain the same + # format expected by CR_FIND. + + if (npasses > 1) { + if (out != in) + call imunmap (out) + call imunmap (in) + im = immap (Memc[output], READ_WRITE, 0) + in = im + out = im + call imseti (in, IM_TYBNDRY, BT_NEAREST) + call imseti (in, IM_NBNDRYPIX, szhwin) + + for (i=2; i<=npasses; i=i+1) { + # Loop through each cosmic ray in the previous pass. + ncr = CR_NCR(cr) + do j = ncrlast+1, ncr { + flag = Memi[CR_FLAG(cr)+j-1] + if (flag==NO || flag==ALWAYSNO) + next + c = Memr[CR_COL(cr)+j-1] + l = Memr[CR_LINE(cr)+j-1] + c1 = max (1-szhwin, c - (szwin-1)) + c2 = min (nc+szhwin, c + (szwin-1)) + k = c2 - c1 + 1 + l1 = max (1-szhwin, l - (szwin-1)) + l2 = min (nl+szhwin, l + (szwin-1)) + + # Set the line pointers off an image section + # centered on a previously replaced cosmic ray. + + ptr = imgs2r (in, c1, c2, l1, l2) - k + + l1 = max (1, l - szhwin) + l2 = min (nl, l + szhwin) + do l = l1, l2 { + do m = 1, szwin + Memi[data+m-1] = ptr + m * k + ptr = ptr + k + + call cr_find ( cr, threshold, Memi[data], + k, szwin, c1, l, sf1, sf2, + Memr[x], Memr[y], Memr[z], Memr[w]) + } + } + call cr_flags (cr, fluxratio) + + # Replace any new cosmic rays found. + call cr_replace (cr, ncr, in, nreplaced) + ncrlast = ncr + } + } + + # Output header log, log, plot, and bad pixels. + call sprintf (Memc[str], SZ_LINE, + "Threshold=%5.1f, fluxratio=%6.2f, removed=%d") + call pargr (threshold) + call pargr (fluxratio) + call pargi (nreplaced) + call imastr (out, "crcor", Memc[str]) + + call cr_plot (cr, in, fluxratio) + call cr_crmask (cr, Memc[badpix], in) + + call cr_close (cr) + if (out != in) + call imunmap (out) + call imunmap (in) + } then { + # In case of error clean up and go on to the next image. + if (in != NULL) { + if (out != NULL && out != in) + call imunmap (out) + call imunmap (in) + } + if (cr != NULL) + call cr_close (cr) + call erract (EA_WARN) + } + } + + if (interactive) { + call gt_free (gt) + call gclose (gp) + } + call imtclose (list1) + call imtclose (list2) + call clpcls (list3) + call gsfree (sf1) + call gsfree (sf2) + call sfree (sp) +end diff --git a/noao/imred/crutil/src/t_craverage.x b/noao/imred/crutil/src/t_craverage.x new file mode 100644 index 00000000..f7b82113 --- /dev/null +++ b/noao/imred/crutil/src/t_craverage.x @@ -0,0 +1,847 @@ +include <error.h> +include <imhdr.h> +include <mach.h> + +define MAXBUF 500000 # Maximum pixel buffer + +define PLSIG 15.87 # Low percentile +define PHSIG 84.13 # High percentile + + +# T_CRAVERAGE -- Detect, fix, and flag cosmic rays. Also detect objects. +# Deviant pixels relative to a local average with the candidate pixel +# excluded and sigma are detected and replaced by the average value +# and/or written to a cosmic ray mask. Average values above a the median +# of a background annulus are detected as objects and cosmic rays are +# excluded. The object positions may be output in the mask. + +procedure t_craverage () + +int inlist # Input image list +int outlist # Output image list +int crlist # Output mask list +int avglist # Output average list +int siglist # Output sigma list +int crval # Output cosmic ray mask value +int objval # Output object mask value +int navg # Averaging box size +int nrej # Number of high pixels to reject from average +int nbkg # Background width +int nsig # Sigma box size +real lobjsig, hobjsig # Object threshold sigmas +real lcrsig, hcrsig # CR threshold sigmas outside of object +real var0 # Variance coefficient for DN^0 term +real var1 # Variance coefficient for DN^1 term +real var2 # Variance coefficient for DN^2 term +real crgrw # Cosmic ray grow radius +real objgrw # Object grow radius + +int i, nc, nl, nlstep, nbox, l1, l2, l3, l4, nl1, pmmode +pointer sp, input, output, crmask, crmask1, extname, average, sigma +pointer in, out, pm, aim, sim +pointer inbuf, pinbuf, outbuf, pbuf, abuf, sbuf + +real clgetr() +int clgeti(), imtopenp(), imtgetim() +pointer immap(), imgs2s(), imgs2r(), imps2r(), imps2s() +errchk immap, imgs2s, imgs2r, imps2r, imps2s, craverage, crgrow, imgstr + +begin + call smark (sp) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (crmask, SZ_FNAME, TY_CHAR) + call salloc (crmask1, SZ_FNAME, TY_CHAR) + call salloc (average, SZ_FNAME, TY_CHAR) + call salloc (sigma, SZ_FNAME, TY_CHAR) + call salloc (extname, SZ_FNAME, TY_CHAR) + + # Get parameters. + inlist = imtopenp ("input") + outlist = imtopenp ("output") + crlist = imtopenp ("crmask") + avglist = imtopenp ("average") + siglist = imtopenp ("sigma") + crval = clgeti ("crval") + objval = clgeti ("objval") + navg = max (1, clgeti ("navg") / 2) + nrej = min (clgeti ("nrej"), navg-1) + nbkg = clgeti ("nbkg") + nsig = clgeti ("nsig") + lobjsig = clgetr ("lobjsig") + hobjsig = clgetr ("hobjsig") + lcrsig = clgetr ("lcrsig") + hcrsig = clgetr ("hcrsig") + nbox = 2 * (navg + nbkg) + 1 + var0 = clgetr ("var0") + var1 = clgetr ("var1") + var2 = clgetr ("var2") + crgrw = clgetr ("crgrow") + objgrw = clgetr ("objgrow") + + # Do the input images. + Memc[crmask1] = EOS + while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) { + if (imtgetim (outlist, Memc[output], SZ_FNAME) == EOF) + Memc[output] = EOS + if (imtgetim (crlist, Memc[crmask], SZ_FNAME) == EOF) + call strcpy (Memc[crmask1], Memc[crmask], SZ_FNAME) + else if (Memc[crmask] == '!') + call strcpy (Memc[crmask], Memc[crmask1], SZ_FNAME) + if (imtgetim (avglist, Memc[average], SZ_FNAME) == EOF) + Memc[average] = EOS + if (imtgetim (siglist, Memc[sigma], SZ_FNAME) == EOF) + Memc[sigma] = EOS + + # Map the input and output images. + iferr { + in = NULL; out = NULL; pm = NULL; aim = NULL; sim = NULL + inbuf = NULL; pinbuf = NULL; outbuf = NULL; pbuf = NULL; + abuf = NULL; sbuf=NULL + + in = immap (Memc[input], READ_ONLY, 0) + if (Memc[output] != EOS) + out = immap (Memc[output], NEW_COPY, in) + if (Memc[crmask] != EOS) { + if (Memc[crmask] == '!') + call imgstr (in, Memc[crmask+1], Memc[crmask], SZ_FNAME) + pmmode = READ_WRITE + iferr (call imgstr (in, "extname", Memc[extname], SZ_FNAME)) + call strcpy ("pl", Memc[extname], SZ_FNAME) + call xt_maskname (Memc[crmask], Memc[extname], pmmode, + Memc[crmask], SZ_FNAME) + iferr (pm = immap (Memc[crmask], pmmode, 0)) { + pmmode = NEW_COPY + pm = immap (Memc[crmask], pmmode, in) + } + } + if (Memc[average] != EOS) + aim = immap (Memc[average], NEW_COPY, in) + if (Memc[sigma] != EOS) + sim = immap (Memc[sigma], NEW_COPY, in) + + # Go through the input in large blocks of lines. If the + # block is smaller than the whole image overlap the blocks + # so the average only has boundaries at the ends of the image. + # However, the output is done in non-overlapping blocks with + # the pointers are adjusted so that addresses can be in the + # space of the input block. CRAVERAGE does not address + # outside of the output data block. Set the mask values + # based on the distances to the nearest good pixels. + + nc = IM_LEN(in,1) + nl = IM_LEN(in,2) + nlstep = max (1, MAXBUF / nc - nbox) + + do i = 1, nl, nlstep { + l1 = i + l2 = min (nl, i + nlstep - 1) + l3 = max (1, l1 - nbox / 2) + l4 = min (nl, l2 + nbox / 2) + nl1 = l4 - l3 + 1 + inbuf = imgs2r (in, 1, nc, l3, l4) + if (out != NULL) + outbuf = imps2r (out, 1, nc, l1, l2) - (l1 - l3) * nc + if (pm != NULL) { + if (pmmode == READ_WRITE) { + pinbuf = imgs2s (pm, 1, nc, l3, l4) + pbuf = imps2s (pm, 1, nc, l1, l2) + call amovs (Mems[pinbuf+(l1-l3)*nc], + Mems[pbuf], nc*(l2-l1+1)) + pbuf = pbuf - (l1 - l3) * nc + } else { + pinbuf = NULL + pbuf = imps2s (pm, 1, nc, l1, l2) + call aclrs (Mems[pbuf], nc*(l2-l1+1)) + pbuf = pbuf - (l1 - l3) * nc + } + } + if (aim != NULL) + abuf = imps2r (aim, 1, nc, l1, l2) - (l1 - l3) * nc + if (sim != NULL) + sbuf = imps2r (sim, 1, nc, l1, l2) - (l1 - l3) * nc + if (pinbuf == NULL) + call craverage (inbuf, outbuf, pbuf, abuf, sbuf, + nc, nl1, l1-l3+1, l2-l3+1, navg, nrej, nbkg, + var0, var1, var2, nsig, lcrsig, hcrsig, + lobjsig, hobjsig, crval, objval) + else + call craverage1 (inbuf, pinbuf, outbuf, pbuf, abuf, + sbuf, nc, nl1, l1-l3+1, l2-l3+1, navg, nrej, nbkg, + var0, var1, var2, nsig, lcrsig, hcrsig, + lobjsig, hobjsig, crval, objval) + } + + # Grow regions if desired. The routines are nops if the + # grow is zero. + + if (pm != NULL) { + if (pmmode != READ_WRITE) { + call imunmap (pm) + iferr (pm = immap (Memc[crmask], READ_WRITE, 0)) + call error (1, "Can't reopen mask for growing") + } + + if (crval == objval) + call crgrow (pm, max (crgrw, objgrw), crval, crval) + else { + call crgrow (pm, crgrw, crval, crval) + call crgrow (pm, objgrw, objval, objval) + } + } + } then + call erract (EA_WARN) + + if (sim != NULL) + call imunmap (sim) + if (aim != NULL) + call imunmap (aim) + if (pm != NULL) + call imunmap (pm) + if (out != NULL) + call imunmap (out) + call imunmap (in) + } + + call imtclose (inlist) + call imtclose (outlist) + call imtclose (crlist) + call imtclose (avglist) + call imtclose (siglist) + + call sfree (sp) +end + + +# CRAVERAGE -- Detect, replace, and flag cosmic rays. +# A local background is computed using moving box averages to avoid +# contaminating bad pixels. If variance model is given then that is +# used otherwise a local sigma is computed in blocks (it is not a moving box +# for efficiency) by using a percentile point of the sorted pixel values to +# estimate the width of the distribution uncontaminated by bad pixels). Once +# the background and sigma are known deviant pixels are found by using sigma +# threshold factors. + +procedure craverage (in, out, pout, aout, sout, nc, nl, nl1, nl2, + navg, nrej, nbkg, var0, var1, var2, nsig, lcrsig, hcrsig, + lobjsig, hobjsig crval, objval) + +pointer in #I Input data +pointer out #O Output data +pointer pout #O Output mask (0=good, 1=bad) +pointer aout #O Output averages +pointer sout #O Output sigmas +int nc, nl #I Number of columns and lines +int nl1, nl2 #I Lines to compute +int navg #I Averaging box half-size +int nrej #I Number of high pixels to reject from average +int nbkg #I Median background width +real var0 #I Variance coefficient for DN^0 term +real var1 #I Variance coefficient for DN^1 term +real var2 #I Variance coefficient for DN^2 term +int nsig #I Sigma box size +real lcrsig, hcrsig #I Threshold sigmas outside of object +real lobjsig, hobjsig #I Object threshold sigmas +int crval #I CR mask value +int objval #I Object mask value + +int i, j, c, c1, c2, c3, c4, l, l1, l2, l3, l4, n1, n2 +int navg2, nbkg2, nsig2, plsig, phsig +real data, avg, bkg, sigma, losig, hosig +real low, high, cravg(), amedr() +pointer stack, avgs, bkgs, sigs, work1, work2 +pointer ptr1, ptr2, ip, op, pp, ap, sp + +begin + navg2 = (2 * navg + 1) ** 2 + nbkg2 = (2 * (navg + nbkg) + 1) ** 2 - navg2 + nsig2 = nsig * nsig + + call smark (stack) + call salloc (avgs, nc, TY_REAL) + call salloc (bkgs, nc, TY_REAL) + call salloc (sigs, nc, TY_REAL) + call salloc (work1, navg2, TY_REAL) + call salloc (work2, max (nsig2, nbkg2), TY_REAL) + + if (var0 != 0. && var1 == 0. && var2 ==0.) + call amovkr (sqrt(var0), Memr[sigs], nc) + + avgs = avgs - 1 + sigs = sigs - 1 + bkgs = bkgs - 1 + + plsig = nint (PLSIG*nsig2/100.-1) + phsig = nint (PHSIG*nsig2/100.-1) + losig = lobjsig / sqrt (real(navg2-1)) + hosig = hobjsig / sqrt (real(navg2-1)) + + do l = nl1, nl2 { + # Compute statistics. + l1 = max (1, l-navg-nbkg) + l2 = max (1, l-navg) + l3 = min (nl, l+navg) + l4 = min (nl, l+navg+nbkg) + ap = aout + (l - 1) * nc + do c = 1, nc { + c1 = max (1, c-navg-nbkg) + c2 = max (1, c-navg) + c3 = min (nc, c+navg) + c4 = min (nc, c+navg+nbkg) + ptr1 = work1 + ptr2 = work2 + n1 = 0 + n2 = 0 + do j = l1, l2-1 { + ip = in + (j - 1) * nc + c1 - 1 + do i = c1, c4 { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + ip = ip + 1 + } + } + do j = l2, l3 { + ip = in + (j - 1) * nc + c1 - 1 + do i = c1, c2-1 { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + ip = ip + 1 + } + do i = c2, c3 { + if (j != l || i != c) { + Memr[ptr1] = Memr[ip] + n1 = n1 + 1 + ptr1 = ptr1 + 1 + } + ip = ip + 1 + } + do i = c3+1, c4 { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + ip = ip + 1 + } + } + do j = l3+1, l4 { + ip = in + (j - 1) * nc + c1 - 1 + do i = c1, c4 { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + ip = ip + 1 + } + } + avg = cravg (Memr[work1], n1, nrej) + bkg = amedr (Memr[work2], n2) + Memr[bkgs+c] = bkg + Memr[avgs+c] = avg + if (aout != NULL) { + Memr[ap] = avg - bkg + ap = ap + 1 + } + } + + # Compute sigmas and output if desired. + if (var0 != 0. || var1 != 0. || var2 != 0.) { + if (var1 != 0.) { + if (var2 != 0.) { + do c = 1, nc { + data = max (0., Memr[avgs+c]) + Memr[sigs+c] = sqrt (var0+var1*data+var2*data**2) + } + } else { + do c = 1, nc { + data = max (0., Memr[avgs+c]) + Memr[sigs+c] = sqrt (var0 + var1 * data) + } + } + } else if (var2 != 0.) { + do c = 1, nc { + data = max (0., Memr[avgs+c]) + Memr[sigs+c] = sqrt (var0 + var2 * data**2) + } + } + } else { + # Compute sigmas from percentiles. This is done in blocks. + if (mod (l-nl1, nsig) == 0 && l<nl-nsig+1) { + do c = 1, nc-nsig+1, nsig { + ptr2 = work2 + n2 = 0 + do j = l, l+nsig-1 { + ip = in + (j - 1) * nc + c - 1 + do i = 1, nsig { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + ip = ip + 1 + } + } + call asrtr (Memr[work2], Memr[work2], n2) + sigma = (Memr[work2+phsig] - Memr[work2+plsig]) / 2. + call amovkr (sigma, Memr[sigs+c], nsig) + } + call amovkr (sigma, Memr[sigs+c], nc-c+1) + } + } + if (sout != NULL) { + sp = sout + (l - 1) * nc + do c = 1, nc { + Memr[sp] = Memr[sigs+c] + sp = sp + 1 + } + } + + # Detect, fix, and flag cosmic rays. + if (pout == NULL && out == NULL) + ; + else if (pout == NULL) { + ip = in + (l - 1) * nc + op = out + (l - 1) * nc + do c = 1, nc { + data = Memr[ip] + avg = Memr[avgs+c] + bkg = Memr[bkgs+c] + sigma = Memr[sigs+c] + low = bkg - losig * sigma + high = bkg + hosig * sigma + if (avg < low || avg > high) { + Memr[op] = data + } else { + low = avg - lcrsig * sigma + high = avg + hcrsig * sigma + if (data < low || data > high) + Memr[op] = avg + else + Memr[op] = data + } + ip = ip + 1 + op = op + 1 + } + } else if (out == NULL) { + ip = in + (l - 1) * nc + pp = pout + (l - 1) * nc + do c = 1, nc { + data = Memr[ip] + avg = Memr[avgs+c] + bkg = Memr[bkgs+c] + sigma = Memr[sigs+c] + low = bkg - losig * sigma + high = bkg + hosig * sigma + if (avg < low || avg > high) + Mems[pp] = objval + else { + low = avg - lcrsig * sigma + high = avg + hcrsig * sigma + if (data < low || data > high) + Mems[pp] = crval + } + ip = ip + 1 + pp = pp + 1 + } + } else { + ip = in + (l - 1) * nc + op = out + (l - 1) * nc + pp = pout + (l - 1) * nc + do c = 1, nc { + data = Memr[ip] + avg = Memr[avgs+c] + bkg = Memr[bkgs+c] + sigma = Memr[sigs+c] + low = bkg - losig * sigma + high = bkg + hosig * sigma + if (avg < low || avg > high) { + Memr[op] = data + Mems[pp] = objval + } else { + low = avg - lcrsig * sigma + high = avg + hcrsig * sigma + if (data < low || data > high) { + Memr[op] = avg + Mems[pp] = crval + } else + Memr[op] = data + } + ip = ip + 1 + op = op + 1 + pp = pp + 1 + } + } + } + + call sfree (stack) +end + + +# CRAVERAGE1 -- Detect, replace, and flag cosmic rays checking input mask. +# A local background is computed using moving box averages to avoid +# contaminating bad pixels. If variance model is given then that is +# used otherwise a local sigma is computed in blocks (it is not a moving box +# for efficiency) by using a percentile point of the sorted pixel values to +# estimate the width of the distribution uncontaminated by bad pixels). Once +# the background and sigma are known deviant pixels are found by using sigma +# threshold factors. + +procedure craverage1 (in, pin, out, pout, aout, sout, nc, nl, nl1, nl2, + navg, nrej, nbkg, var0, var1, var2, nsig, lcrsig, hcrsig, + lobjsig, hobjsig crval, objval) + +pointer in #I Input data +pointer pin #I Pixel mask data +pointer out #O Output data +pointer pout #O Output mask (0=good, 1=bad) +pointer aout #O Output averages +pointer sout #O Output sigmas +int nc, nl #I Number of columns and lines +int nl1, nl2 #I Lines to compute +int navg #I Averaging box half-size +int nrej #I Number of high pixels to reject from average +int nbkg #I Median background width +real var0 #I Variance coefficient for DN^0 term +real var1 #I Variance coefficient for DN^1 term +real var2 #I Variance coefficient for DN^2 term +int nsig #I Sigma box size +real lcrsig, hcrsig #I Threshold sigmas outside of object +real lobjsig, hobjsig #I Object threshold sigmas +int crval #I CR mask value +int objval #I Object mask value + +int i, j, c, c1, c2, c3, c4, l, l1, l2, l3, l4, n1, n2 +int navg2, nbkg2, nsig2, plsig, phsig +real data, avg, bkg, sigma, losig, hosig +real low, high, cravg(), amedr() +pointer stack, avgs, bkgs, sigs, work1, work2 +pointer ptr1, ptr2, ip, mp, op, pp, ap, sp + +begin + navg2 = (2 * navg + 1) ** 2 + nbkg2 = (2 * (navg + nbkg) + 1) ** 2 - navg2 + nsig2 = nsig * nsig + + call smark (stack) + call salloc (avgs, nc, TY_REAL) + call salloc (bkgs, nc, TY_REAL) + call salloc (sigs, nc, TY_REAL) + call salloc (work1, navg2, TY_REAL) + call salloc (work2, max (nsig2, nbkg2), TY_REAL) + + if (var0 != 0. && var1 == 0. && var2 ==0.) + call amovkr (sqrt(var0), Memr[sigs], nc) + + avgs = avgs - 1 + sigs = sigs - 1 + bkgs = bkgs - 1 + + losig = lobjsig / sqrt (real(navg2-1)) + hosig = hobjsig / sqrt (real(navg2-1)) + + do l = nl1, nl2 { + # Compute statistics. + l1 = max (1, l-navg-nbkg) + l2 = max (1, l-navg) + l3 = min (nl, l+navg) + l4 = min (nl, l+navg+nbkg) + ap = aout + (l - 1) * nc + do c = 1, nc { + c1 = max (1, c-navg-nbkg) + c2 = max (1, c-navg) + c3 = min (nc, c+navg) + c4 = min (nc, c+navg+nbkg) + ptr1 = work1 + ptr2 = work2 + n1 = 0 + n2 = 0 + do j = l1, l2-1 { + ip = in + (j - 1) * nc + c1 - 1 + mp = pin + (j - 1) * nc + c1 - 1 + do i = c1, c4 { + if (Mems[mp] == 0) { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + } + ip = ip + 1 + mp = mp + 1 + } + } + do j = l2, l3 { + ip = in + (j - 1) * nc + c1 - 1 + mp = pin + (j - 1) * nc + c1 - 1 + do i = c1, c2-1 { + if (Mems[mp] == 0) { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + } + ip = ip + 1 + mp = mp + 1 + } + do i = c2, c3 { + if ((j != l || i != c) && Mems[mp] == 0) { + Memr[ptr1] = Memr[ip] + n1 = n1 + 1 + ptr1 = ptr1 + 1 + } + ip = ip + 1 + mp = mp + 1 + } + do i = c3+1, c4 { + if (Mems[mp] == 0) { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + } + ip = ip + 1 + mp = mp + 1 + } + } + do j = l3+1, l4 { + ip = in + (j - 1) * nc + c1 - 1 + mp = pin + (j - 1) * nc + c1 - 1 + do i = c1, c4 { + if (Mems[mp] == 0) { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + } + ip = ip + 1 + } + } + if (n1 > 0) + avg = cravg (Memr[work1], n1, nrej) + else + avg = INDEFR + if (n2 > 0) + bkg = amedr (Memr[work2], n2) + else + bkg = INDEFR + Memr[bkgs+c] = bkg + Memr[avgs+c] = avg + if (aout != NULL) { + if (IS_INDEFR(avg) || IS_INDEFR(bkg)) + Memr[ap] = 0. + else + Memr[ap] = avg - bkg + ap = ap + 1 + } + } + + # Compute sigmas and output if desired. + if (var0 != 0. || var1 != 0. || var2 != 0.) { + if (var1 != 0.) { + if (var2 != 0.) { + do c = 1, nc { + data = max (0., Memr[avgs+c]) + Memr[sigs+c] = sqrt (var0+var1*data+var2*data**2) + } + } else { + do c = 1, nc { + data = max (0., Memr[avgs+c]) + Memr[sigs+c] = sqrt (var0 + var1 * data) + } + } + } else if (var2 != 0.) { + do c = 1, nc { + data = max (0., Memr[avgs+c]) + Memr[sigs+c] = sqrt (var0 + var2 * data**2) + } + } + } else { + # Compute sigmas from percentiles. This is done in blocks. + if (mod (l-nl1, nsig) == 0 && l<nl-nsig+1) { + do c = 1, nc-nsig+1, nsig { + ptr2 = work2 + n2 = 0 + do j = l, l+nsig-1 { + ip = in + (j - 1) * nc + c - 1 + mp = pin + (j - 1) * nc + c - 1 + do i = 1, nsig { + if (Mems[mp] == 0) { + Memr[ptr2] = Memr[ip] + n2 = n2 + 1 + ptr2 = ptr2 + 1 + } + ip = ip + 1 + mp = mp + 1 + } + } + if (n2 > 10) { + call asrtr (Memr[work2], Memr[work2], n2) + plsig = nint (PLSIG*n2/100.-1) + phsig = nint (PHSIG*n2/100.-1) + sigma = (Memr[work2+phsig]-Memr[work2+plsig])/2. + } else + sigma = INDEFR + call amovkr (sigma, Memr[sigs+c], nsig) + } + call amovkr (sigma, Memr[sigs+c], nc-c+1) + } + } + if (sout != NULL) { + sp = sout + (l - 1) * nc + do c = 1, nc { + sigma = Memr[sigs+c] + if (IS_INDEFR(sigma)) + Memr[sp] = 0. + else + Memr[sp] = sigma + sp = sp + 1 + } + } + + # Detect, fix, and flag cosmic rays. + if (pout == NULL && out == NULL) + ; + if (pout == NULL) { + ip = in + (l - 1) * nc + mp = pin + (l - 1) * nc + op = out + (l - 1) * nc + do c = 1, nc { + data = Memr[ip] + avg = Memr[avgs+c] + bkg = Memr[bkgs+c] + sigma = Memr[sigs+c] + if (!(Mems[mp] != 0 || IS_INDEFR(avg) || + IS_INDEFR(bkg) || IS_INDEFR(sigma))) { + low = bkg - losig * sigma + high = bkg + hosig * sigma + if (avg < low || avg > high) { + Memr[op] = data + } else { + low = avg - lcrsig * sigma + high = avg + hcrsig * sigma + if (data < low || data > high) + Memr[op] = avg + else + Memr[op] = data + } + } else + Memr[op] = data + ip = ip + 1 + mp = mp + 1 + op = op + 1 + } + } else if (out == NULL) { + ip = in + (l - 1) * nc + mp = pin + (l - 1) * nc + pp = pout + (l - 1) * nc + do c = 1, nc { + data = Memr[ip] + avg = Memr[avgs+c] + bkg = Memr[bkgs+c] + sigma = Memr[sigs+c] + if (!(Mems[mp] != 0 || IS_INDEFR(avg) || + IS_INDEFR(bkg) || IS_INDEFR(sigma))) { + low = bkg - losig * sigma + high = bkg + hosig * sigma + if (avg < low || avg > high) + Mems[pp] = objval + else { + low = avg - lcrsig * sigma + high = avg + hcrsig * sigma + if (data < low || data > high) + Mems[pp] = crval + } + } + ip = ip + 1 + mp = mp + 1 + pp = pp + 1 + } + } else { + ip = in + (l - 1) * nc + mp = pin + (l - 1) * nc + op = out + (l - 1) * nc + pp = pout + (l - 1) * nc + do c = 1, nc { + data = Memr[ip] + avg = Memr[avgs+c] + bkg = Memr[bkgs+c] + sigma = Memr[sigs+c] + if (!(Mems[mp] != 0 || IS_INDEFR(avg) || + IS_INDEFR(bkg) || IS_INDEFR(sigma))) { + low = bkg - losig * sigma + high = bkg + hosig * sigma + if (avg < low || avg > high) { + Memr[op] = data + Mems[pp] = objval + } else { + low = avg - lcrsig * sigma + high = avg + hcrsig * sigma + if (data < low || data > high) { + Memr[op] = avg + Mems[pp] = crval + } else + Memr[op] = data + } + } else + Memr[op] = data + ip = ip + 1 + mp = mp + 1 + op = op + 1 + pp = pp + 1 + } + } + } + + call sfree (stack) +end + + +# CRAVG -- Compute average with the highest nrej points excluded. +# When nrej is greater than 2 the data array will be returned sorted. + +real procedure cravg (data, npts, nrej) + +real data[npts] #I Input data (will be sorted if nrej>2) +int npts #I Number of data points +int nrej #I Number of data points to reject + +int i +real sum, max1, max2, val + +begin + if (npts <= nrej) + return (INDEFR) + + switch (nrej) { + case 0: + sum = 0. + do i = 1, npts + sum = sum + data[i] + case 1: + sum = 0. + max1 = data[1] + do i = 2, npts { + val = data[i] + if (val > max1) { + sum = sum + max1 + max1 = val + } else + sum = sum + val + } + case 2: + sum = 0. + max1 = min (data[1], data[2]) + max2 = max (data[1], data[2]) + do i = 3, npts { + val = data[i] + if (val > max1) { + sum = sum + max1 + if (val > max2) { + max1 = max2 + max2 = val + } else + max1 = val + } else + sum = sum + val + } + default: + call asrtr (data, data, npts) + sum = 0. + do i = 1, npts-nrej + sum = sum + data[i] + } + + return (sum / (npts - nrej)) +end diff --git a/noao/imred/crutil/src/t_crgrow.x b/noao/imred/crutil/src/t_crgrow.x new file mode 100644 index 00000000..7316cc76 --- /dev/null +++ b/noao/imred/crutil/src/t_crgrow.x @@ -0,0 +1,182 @@ +include <error.h> +include <imhdr.h> + +# T_CRGROW -- Grow cosmic ray mask identifications. + +procedure t_crgrow () + +int input # Input masks +int output # Output masks +real radius # Radius +int inval # Input mask value to grow +int outval # Output grown mask value + +pointer sp, inmask, outmask, temp1, temp2 +pointer im, ptr + +int imtopenp(), imtlen(), imtgetim() +bool strne() +int clgeti() +real clgetr() +pointer immap() +errchk immap, crgrow + +begin + call smark (sp) + call salloc (inmask, SZ_FNAME, TY_CHAR) + call salloc (outmask, SZ_FNAME, TY_CHAR) + call salloc (temp1, SZ_FNAME, TY_CHAR) + call salloc (temp2, SZ_FNAME, TY_CHAR) + + # Task parameters. + input = imtopenp ("input") + output = imtopenp ("output") + radius = max (0., clgetr ("radius")) + inval = clgeti ("inval") + outval = clgeti ("outval") + + if (imtlen (output) != imtlen (input)) + call error (1, "Input and output lists do not match") + + # Grow the cosmic ray masks. + while (imtgetim (input, Memc[inmask], SZ_FNAME) != EOF) { + call strcpy (Memc[inmask], Memc[outmask], SZ_FNAME) + if (imtgetim (output, Memc[outmask], SZ_FNAME) == EOF) + call error (1, "Output list ended prematurely") + if (strne (Memc[inmask], Memc[outmask])) { + call imgcluster (Memc[inmask], Memc[temp1], SZ_FNAME) + call imgcluster (Memc[outmask], Memc[temp2], SZ_FNAME) + iferr (call imcopy (Memc[temp1], Memc[temp2])) { + call erract (EA_WARN) + next + } + im = immap (Memc[inmask], READ_ONLY, TY_CHAR) + iferr (call imgstr (im, "extname", Memc[temp1], SZ_FNAME)) + call strcpy ("pl", Memc[temp1], SZ_FNAME) + call imunmap (im) + } + call xt_maskname (Memc[outmask], Memc[temp1], 0, Memc[outmask], + SZ_FNAME) + + if (radius < 1.) + next + + iferr { + im = NULL + ptr = immap (Memc[outmask], READ_WRITE, 0); im = ptr + call crgrow (im, radius, inval, outval) + } then { + call erract (EA_WARN) + if (strne (Memc[inmask], Memc[outmask])) { + if (im != NULL) { + call imunmap (im) + iferr (call imdelete (Memc[outmask])) + call erract (EA_WARN) + } + } + } + + if (im != NULL) + call imunmap (im) + } + + call imtclose (output) + call imtclose (input) + call sfree (sp) +end + + +# CRGROW -- Grow cosmic rays. + +procedure crgrow (im, grow, inval, outval) + +pointer im # Mask pointer (Read/Write) +real grow # Radius (pixels) +int inval # Input mask value for pixels to grow +int outval # Output mask value for grown pixels + +int i, j, k, l, nc, nl, ngrow, nbufs, val1, val2 +long v1[2], v2[2] +real grow2, y2 +pointer , buf, buf1, buf2, ptr + +int imgnli(), impnli() +errchk calloc, imgnli, impnli + +begin + if (grow < 1. || inval == 0) + return + + grow2 = grow * grow + ngrow = int (grow) + buf = NULL + + iferr { + if (IM_NDIM(im) > 2) + call error (1, + "Only one or two dimensional masks are allowed") + + nc = IM_LEN(im, 1) + if (IM_NDIM(im) > 1) + nl = IM_LEN(im,2) + else + nl = 1 + + # Initialize buffering. + nbufs = min (1 + 2 * ngrow, nl) + call calloc (buf, nc*nbufs, TY_INT) + + call amovkl (long(1), v1, IM_NDIM(im)) + call amovkl (long(1), v2, IM_NDIM(im)) + while (imgnli (im, buf1, v1) != EOF) { + j = v1[2] - 1 + buf2 = buf + mod (j, nbufs) * nc + do i = 1, nc { + val1 = Memi[buf1] + val2 = Memi[buf2] + if ((IS_INDEFI(inval) && val1 != 0) || val1 == inval) { + do k = max(1,j-ngrow), min (nl,j+ngrow) { + ptr = buf + mod (k, nbufs) * nc - 1 + y2 = (k - j) ** 2 + do l = max(1,i-ngrow), min (nc,i+ngrow) { + if ((l-i)**2 + y2 > grow2) + next + Memi[ptr+l] = -val1 + } + } + } else { + if (val2 >= 0) + Memi[buf2] = val1 + } + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + + if (j > ngrow) { + while (impnli (im, buf2, v2) != EOF) { + k = v2[2] - 1 + buf1 = buf + mod (k, nbufs) * nc + do i = 1, nc { + val1 = Memi[buf1] + if (val1 < 0) { + if (IS_INDEFI(outval)) + Memi[buf2] = -val1 + else + Memi[buf2] = outval + } else + Memi[buf2] = val1 + Memi[buf1] = 0. + buf1 = buf1 + 1 + buf2 = buf2 + 1 + } + if (j != nl) + break + } + } + } + } then + call erract (EA_ERROR) + + if (buf != NULL) + call mfree (buf, TY_INT) +end diff --git a/noao/imred/crutil/src/t_crmedian.x b/noao/imred/crutil/src/t_crmedian.x new file mode 100644 index 00000000..6c7d0fba --- /dev/null +++ b/noao/imred/crutil/src/t_crmedian.x @@ -0,0 +1,417 @@ +include <imhdr.h> +include <mach.h> + +define MAXBUF 500000 # Maximum pixel buffer + +define PLSIG 15.87 # Low percentile +define PHSIG 84.13 # High percentile + + +# T_CRMEDIAN -- Detect, fix, and flag cosmic rays. +# Deviant pixels relative to a local median and sigma are detected and +# replaced by the median value and/or written to a cosmic ray mask. + +procedure t_crmedian () + +pointer input # Input image +pointer output # Output image +pointer crmask # Output mask +pointer median # Output median +pointer sigma # Output sigma +pointer residual # Output residual +real var0 # Variance coefficient for DN^0 term +real var1 # Variance coefficient for DN^1 term +real var2 # Variance coefficient for DN^2 term +real lsig, hsig # Threshold sigmas +int ncmed, nlmed # Median box size +int ncsig, nlsig # Sigma box size + +int i, nc, nl, nlstep, l1, l2, l3, l4, nl1 +pointer sp, extname, in, out, pm, mim, sim, rim +pointer inbuf, outbuf, pmbuf, mbuf, sbuf, rbuf +real clgetr() +int clgeti(), nowhite() +pointer immap(), imgs2r(), imps2r(), imps2s() +errchk immap, imgs2r, imps2r, imps2s, crmedian, imgstr + +begin + call smark (sp) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (output, SZ_FNAME, TY_CHAR) + call salloc (crmask, SZ_FNAME, TY_CHAR) + call salloc (residual, SZ_FNAME, TY_CHAR) + call salloc (median, SZ_FNAME, TY_CHAR) + call salloc (sigma, SZ_FNAME, TY_CHAR) + call salloc (extname, SZ_FNAME, TY_CHAR) + + # Get parameters. + call clgstr ("input", Memc[input], SZ_FNAME) + call clgstr ("output", Memc[output], SZ_FNAME) + call clgstr ("crmask", Memc[crmask], SZ_FNAME) + call clgstr ("median", Memc[median], SZ_FNAME) + call clgstr ("sigma", Memc[sigma], SZ_FNAME) + call clgstr ("residual", Memc[residual], SZ_FNAME) + var0 = clgetr ("var0") + var1 = clgetr ("var1") + var2 = clgetr ("var2") + lsig = clgetr ("lsigma") + hsig = clgetr ("hsigma") + ncmed = clgeti ("ncmed") + nlmed = clgeti ("nlmed") + ncsig = clgeti ("ncsig") + nlsig = clgeti ("nlsig") + + # Map the input and output images. + in = NULL; out = NULL; pm = NULL; mim = NULL; sim = NULL; rim = NULL + inbuf = NULL; outbuf = NULL; pmbuf = NULL + mbuf = NULL; sbuf=NULL; rbuf = NULL + in = immap (Memc[input], READ_ONLY, 0) + if (nowhite (Memc[output], Memc[output], SZ_FNAME) > 0) + out = immap (Memc[output], NEW_COPY, in) + if (nowhite (Memc[crmask], Memc[crmask], SZ_FNAME) > 0) { + if (Memc[crmask] == '!') + call imgstr (in, Memc[crmask+1], Memc[crmask], SZ_FNAME) + iferr (call imgstr (in, "extname", Memc[extname], SZ_FNAME)) + call strcpy ("pl", Memc[extname], SZ_FNAME) + call xt_maskname (Memc[crmask], Memc[extname], 0, Memc[crmask], + SZ_FNAME) + pm = immap (Memc[crmask], NEW_COPY, in) + } + if (nowhite (Memc[median], Memc[median], SZ_FNAME) > 0) + mim = immap (Memc[median], NEW_COPY, in) + if (nowhite (Memc[sigma], Memc[sigma], SZ_FNAME) > 0) + sim = immap (Memc[sigma], NEW_COPY, in) + if (nowhite (Memc[residual], Memc[residual], SZ_FNAME) > 0) + rim = immap (Memc[residual], NEW_COPY, in) + + # Go through the input in large blocks of lines. If the + # block is smaller than the whole image overlap the blocks + # so the median only has boundaries at the ends of the image. + # However, the output is done in non-overlapping blocks with + # the pointers are adjusted so that addresses can be in the space + # of the input block. CRMEDIAN does not address outside of + # the output data block. Set the mask values based on the + # distances to the nearest good pixels. + + nc = IM_LEN(in,1) + nl = IM_LEN(in,2) + nlstep = max (1, MAXBUF / nc - nlmed) + + do i = 1, nl, nlstep { + l1 = i + l2 = min (nl, i + nlstep - 1) + l3 = max (1, l1 - nlmed / 2) + l4 = min (nl, l2 + nlmed / 2) + nl1 = l4 - l3 + 1 + inbuf = imgs2r (in, 1, nc, l3, l4) + if (out != NULL) + outbuf = imps2r (out, 1, nc, l1, l2) - (l1 - l3) * nc + if (pm != NULL) + pmbuf = imps2s (pm, 1, nc, l1, l2) - (l1 - l3) * nc + if (mim != NULL) + mbuf = imps2r (mim, 1, nc, l1, l2) - (l1 - l3) * nc + if (sim != NULL) + sbuf = imps2r (sim, 1, nc, l1, l2) - (l1 - l3) * nc + if (rim != NULL) + rbuf = imps2r (rim, 1, nc, l1, l2) - (l1 - l3) * nc + call crmedian (inbuf, outbuf, pmbuf, mbuf, sbuf, rbuf, + nc, nl1, l1-l3+1, l2-l3+1, ncmed, nlmed, var0, var1, var2, + ncsig, nlsig, lsig, hsig) + } + + if (rim != NULL) + call imunmap (rim) + if (sim != NULL) + call imunmap (sim) + if (mim != NULL) + call imunmap (mim) + if (pm != NULL) + call imunmap (pm) + if (out != NULL) + call imunmap (out) + call imunmap (in) + call sfree (sp) +end + + +# CRMEDIAN -- Detect, replace, and flag cosmic rays. +# A local background is computed using moving box medians to avoid +# contaminating bad pixels. If variance model is given then that is +# used otherwise a local sigma is computed in blocks (it is not a moving box +# for efficiency) by using a percentile point of the sorted pixel values to +# estimate the width of the distribution uncontaminated by bad pixels). Once +# the background and sigma are known deviant pixels are found by using sigma +# threshold factors. + +procedure crmedian (in, out, pout, mout, sout, rout, nc, nl, nl1, nl2, + ncmed, nlmed, var0, var1, var2, ncsig, nlsig, lsig, hsig) + +pointer in #I Input data +pointer out #O Output data +pointer pout #O Output mask (0=good, 1=bad) +pointer mout #O Output medians +pointer sout #O Output sigmas +pointer rout #O Output residuals +int nc, nl #I Number of columns and lines +int nl1, nl2 #I Lines to compute +int ncmed, nlmed #I Median box size +real var0 # Variance coefficient for DN^0 term +real var1 # Variance coefficient for DN^1 term +real var2 # Variance coefficient for DN^2 term +int ncsig, nlsig #I Sigma box size +real lsig, hsig #I Threshold sigmas + +int i, j, k, l, m, plsig, phsig +real data, med, sigma, low, high, amedr() +pointer stack, meds, sigs, work, ptr, ip, op, pp, mp, sp, rp + +begin + call smark (stack) + call salloc (meds, nc, TY_REAL) + call salloc (sigs, nc, TY_REAL) + call salloc (work, max (ncsig*nlsig, ncmed*nlmed), TY_REAL) + + if (var0 != 0. && var1 == 0. && var2 ==0.) + call amovkr (sqrt(var0), Memr[sigs], nc) + + meds = meds - 1 + sigs = sigs - 1 + + i = ncsig * nlsig + plsig = nint (PLSIG*i/100.-1) + phsig = nint (PHSIG*i/100.-1) + + do i = nl1, nl2 { + + # Compute median and output if desired. This is a moving median. + l = min (nl, i+nlmed/2) + l = max (1, l-nlmed+1) + mp = mout + (i - 1) * nc + do j = 1, nc { + k = min (nc, j+ncmed/2) + k = max (1, k-ncmed+1) + ptr = work + ip = in + (l - 1) * nc + k - 1 + do m = l, l+nlmed-1 { + call amovr (Memr[ip], Memr[ptr], ncmed) + ip = ip + nc + ptr = ptr + ncmed + } + med = amedr (Memr[work], ncmed * nlmed) + Memr[meds+j] = med + if (mout != NULL) { + Memr[mp] = med + mp = mp + 1 + } + } + + # Compute sigmas and output if desired. + if (var0 != 0. || var1 != 0. || var2 != 0.) { + if (var1 != 0.) { + if (var2 != 0.) { + do j = 1, nc { + data = max (0., Memr[meds+j]) + Memr[sigs+j] = sqrt (var0 + var1*data + var2*data**2) + } + } else { + do j = 1, nc { + data = max (0., Memr[meds+j]) + Memr[sigs+j] = sqrt (var0 + var1 * data) + } + } + } else if (var2 != 0.) { + do j = 1, nc { + data = max (0., Memr[meds+j]) + Memr[sigs+j] = sqrt (var0 + var2 * data**2) + } + } + } else { + # Compute sigmas from percentiles. This is done in blocks. + if (mod (i-nl1, nlsig) == 0 && i<nl-nlsig+1) { + do j = 1, nc-ncsig+1, ncsig { + ptr = work + ip = in + (i - 1) * nc + j - 1 + do k = i, i+nlsig-1 { + call amovr (Memr[ip], Memr[ptr], ncsig) + ip = ip + nc + ptr = ptr + ncsig + } + call asrtr (Memr[work], Memr[work], ncsig*nlsig) + sigma = (Memr[work+phsig] - Memr[work+plsig]) / 2. + call amovkr (sigma, Memr[sigs+j], ncsig) + } + call amovkr (sigma, Memr[sigs+j], nc-j+1) + } + } + if (sout != NULL) { + sp = sout + (i - 1) * nc + do j = 1, nc { + Memr[sp] = Memr[sigs+j] + sp = sp + 1 + } + } + + # Detect, fix, and flag cosmic rays. + if (rout == NULL) { + if (pout == NULL) { + ip = in + (i - 1) * nc + op = out + (i - 1) * nc + do j = 1, nc { + data = Memr[ip] + med = Memr[meds+j] + sigma = Memr[sigs+j] + low = med - lsig * sigma + high = med + hsig * sigma + if (data < low || data > high) + Memr[op] = med + else + Memr[op] = data + ip = ip + 1 + op = op + 1 + } + } else if (out == NULL) { + ip = in + (i - 1) * nc + pp = pout + (i - 1) * nc + do j = 1, nc { + data = Memr[ip] + med = Memr[meds+j] + sigma = Memr[sigs+j] + low = med - lsig * sigma + high = med + hsig * sigma + if (data < low || data > high) + Mems[pp] = 1 + else + Mems[pp] = 0 + ip = ip + 1 + pp = pp + 1 + } + } else { + ip = in + (i - 1) * nc + op = out + (i - 1) * nc + pp = pout + (i - 1) * nc + do j = 1, nc { + data = Memr[ip] + med = Memr[meds+j] + sigma = Memr[sigs+j] + low = med - lsig * sigma + high = med + hsig * sigma + if (data < low || data > high) { + Memr[op] = med + Mems[pp] = 1 + } else { + Memr[op] = data + Mems[pp] = 0 + } + ip = ip + 1 + op = op + 1 + pp = pp + 1 + } + } + } else { + if (pout == NULL && out == NULL) { + ip = in + (i - 1) * nc + rp = rout + (i - 1) * nc + do j = 1, nc { + data = Memr[ip] + med = Memr[meds+j] + sigma = Memr[sigs+j] + if (sigma > 0.) + Memr[rp] = (data - med) / sigma + else { + if ((data - med) < 0.) + Memr[rp] = -MAX_REAL + else + Memr[rp] = MAX_REAL + } + ip = ip + 1 + rp = rp + 1 + } + } else if (pout == NULL) { + ip = in + (i - 1) * nc + op = out + (i - 1) * nc + rp = rout + (i - 1) * nc + do j = 1, nc { + data = Memr[ip] + med = Memr[meds+j] + sigma = Memr[sigs+j] + low = med - lsig * sigma + high = med + hsig * sigma + if (data < low || data > high) + Memr[op] = med + else + Memr[op] = data + if (sigma > 0.) + Memr[rp] = (data - med) / sigma + else { + if ((data - med) < 0.) + Memr[rp] = -MAX_REAL + else + Memr[rp] = MAX_REAL + } + ip = ip + 1 + op = op + 1 + rp = rp + 1 + } + } else if (out == NULL) { + ip = in + (i - 1) * nc + pp = pout + (i - 1) * nc + rp = rout + (i - 1) * nc + do j = 1, nc { + data = Memr[ip] + med = Memr[meds+j] + sigma = Memr[sigs+j] + low = med - lsig * sigma + high = med + hsig * sigma + if (data < low || data > high) + Mems[pp] = 1 + else + Mems[pp] = 0 + if (sigma > 0.) + Memr[rp] = (data - med) / sigma + else { + if ((data - med) < 0.) + Memr[rp] = -MAX_REAL + else + Memr[rp] = MAX_REAL + } + ip = ip + 1 + pp = pp + 1 + rp = rp + 1 + } + } else { + ip = in + (i - 1) * nc + op = out + (i - 1) * nc + pp = pout + (i - 1) * nc + rp = rout + (i - 1) * nc + do j = 1, nc { + data = Memr[ip] + med = Memr[meds+j] + sigma = Memr[sigs+j] + low = med - lsig * sigma + high = med + hsig * sigma + if (data < low || data > high) { + Memr[op] = med + Mems[pp] = 1 + } else { + Memr[op] = data + Mems[pp] = 0 + } + if (sigma > 0.) + Memr[rp] = (data - med) / sigma + else { + if ((data - med) < 0.) + Memr[rp] = -MAX_REAL + else + Memr[rp] = MAX_REAL + } + ip = ip + 1 + op = op + 1 + pp = pp + 1 + rp = rp + 1 + } + } + } + } + + call sfree (stack) +end diff --git a/noao/imred/crutil/src/x_crutil.x b/noao/imred/crutil/src/x_crutil.x new file mode 100644 index 00000000..48408a9f --- /dev/null +++ b/noao/imred/crutil/src/x_crutil.x @@ -0,0 +1,4 @@ +task cosmicrays = t_cosmicrays, + craverage = t_craverage, + crgrow = t_crgrow, + crmedian = t_crmedian diff --git a/noao/imred/crutil/src/xtmaskname.x b/noao/imred/crutil/src/xtmaskname.x new file mode 100644 index 00000000..eb132e12 --- /dev/null +++ b/noao/imred/crutil/src/xtmaskname.x @@ -0,0 +1,125 @@ +# MASKNAME -- Make a mask name. This creates a FITS mask extension if +# possible, otherwise it creates a pixel list file. To create a FITS +# extension the filename must explicitly select the FITS kernel or the +# default image type must be a FITS file. The input and output strings +# may be the same. + +procedure xt_maskname (fname, extname, mode, mname, maxchar) + +char fname[ARB] #I File name +char extname[ARB] #I Default pixel mask extension name +int mode #I Mode +char mname[maxchar] #O Output mask name +int maxchar #I Maximum characters in mask name + +int i, fits +pointer sp, temp + +bool streq() +int strmatch(), stridxs(), strldxs(), strncmp() +int envfind(), access(), imaccess() + +begin + call smark (sp) + call salloc (temp, maxchar, TY_CHAR) + + # Determine whether to use FITS pixel mask extensions. One may set + # fits=NO to force use of pl even when FITS mask extensions are + # supported. + fits = YES + if (fits == YES && envfind ("masktype", Memc[temp], maxchar) > 0) { + if (streq (Memc[temp], "pl")) + fits = NO + } + i = strldxs ("]", fname) + + # Check for explicit .pl extension. + if (strmatch (fname, ".pl$") > 0) + call strcpy (fname, mname, maxchar) + + # Check for explicit mask extension. + else if (strmatch (fname, "type=mask") > 0) + call strcpy (fname, mname, maxchar) + else if (strmatch (fname, "type\\\=mask") > 0) + call strcpy (fname, mname, maxchar) + + # Add mask type. + else if (i > 0) { + if (mode != READ_ONLY) { + call strcpy (fname[i], Memc[temp], maxchar) + call sprintf (mname[i], maxchar-i, ",type=mask%s") + call pargstr (Memc[temp]) + } + i = stridxs ("[", mname) + if (i > 0) { + call strcpy (mname[i+1], Memc[temp], SZ_FNAME) + if (strncmp (mname[i+1], "append", 6)==0 || + strncmp (mname[i+1], "inherit", 7)==0) { + call sprintf (mname[i+1], SZ_FNAME-i+1, "%s,%s") + call pargstr (extname) + call pargstr (Memc[temp]) + } + } + + # Create output from rootname name. + } else if (fits == YES) { + call strcpy (fname, Memc[temp], SZ_FNAME) + if (mode == READ_ONLY) { + call sprintf (mname, maxchar, "%s[%s]") + call pargstr (Memc[temp]) + call pargstr (extname) + } else { + call sprintf (mname, maxchar, "%s[%s,type=mask]") + call pargstr (Memc[temp]) + call pargstr (extname) + } + } else + call strcat (".pl", mname, maxchar) + + # Check if extension name is actually given. + i = stridxs ("[", mname) + if (i > 0) { + call strcpy (mname[i+1], Memc[temp], SZ_FNAME) + if (strncmp (mname[i+1], "append", 6)==0 || + strncmp (mname[i+1], "inherit", 7)==0) { + call sprintf (mname[i+1], SZ_FNAME-i+1, "%s,%s") + call pargstr (extname) + call pargstr (Memc[temp]) + } + } + + # Convert to pl form if required. + i = stridxs ("[", mname) + if (i > 0 && mode == READ_ONLY) + fits = imaccess (mname, mode) + if (fits == NO && i > 0) { + mname[i] = EOS + if (mode == NEW_IMAGE) { + if (access (mname, 0, 0) == NO) { + ifnoerr (call fmkdir (mname)) + mname[i] = '/' + else + mname[i] = '.' + } else + mname[i] = '/' + } else { + if (access (mname, 0, 0) == NO) + mname[i] = '.' + else + mname[i] = '/' + } + + if (strncmp (mname[i+1], "type", 4) == 0 || + strncmp (mname[i+1], "append", 6) == 0 || + strncmp (mname[i+1], "inherit", 7) == 0) { + mname[i+1] = EOS + call strcat (extname, mname, maxchar) + } else { + i = stridxs (",]", mname) + mname[i] = EOS + } + call strcat (".pl", mname, maxchar) + } + + call sfree (sp) +end |