diff options
Diffstat (limited to 'noao/imred/crutil')
38 files changed, 4985 insertions, 0 deletions
diff --git a/noao/imred/crutil/crutil.cl b/noao/imred/crutil/crutil.cl new file mode 100644 index 00000000..89e7ff3e --- /dev/null +++ b/noao/imred/crutil/crutil.cl @@ -0,0 +1,18 @@ +#{ CRUTIL -- Cosmic Ray Utility Package + +package crutil + +reset crusrc = "crutil$src/" + +task crcombine = crusrc$crcombine.cl +task crnebula = crusrc$crnebula.cl +task crfix = crusrc$crfix.cl +task credit = crusrc$credit.cl + +task cosmicrays, + craverage, + crgrow, + crmedian = crusrc$x_crutil.e + + +clbye() diff --git a/noao/imred/crutil/crutil.hd b/noao/imred/crutil/crutil.hd new file mode 100644 index 00000000..0b6ee225 --- /dev/null +++ b/noao/imred/crutil/crutil.hd @@ -0,0 +1,17 @@ +# Help directory for the CRUTIL package. + +$doc = "./doc/" +$crusrc = "./src/" + +revisions sys=crusrc$Revisions + +overview hlp=doc$overview.hlp + +cosmicrays hlp=doc$cosmicrays.hlp, src=crusrc$t_cosmicrays.x +craverage hlp=doc$craverage.hlp, src=crusrc$t_craverage.x +crcombine hlp=doc$crcombine.hlp, src=crusrc$crcombine.cl +credit hlp=doc$credit.hlp, src=crusrc$credit.cl +crfix hlp=doc$crfix.hlp, src=crusrc$crfix.cl +crgrow hlp=doc$crgrow.hlp, src=crusrc$t_crgrow.x +crmedian hlp=doc$crmedian.hlp, src=crusrc$t_crmedian.x +crnebula hlp=doc$crnebula.hlp, src=crusrc$crnebula.cl diff --git a/noao/imred/crutil/crutil.men b/noao/imred/crutil/crutil.men new file mode 100644 index 00000000..f6117e5f --- /dev/null +++ b/noao/imred/crutil/crutil.men @@ -0,0 +1,10 @@ + cosmicrays - Remove cosmic rays using flux ratio algorithm + craverage - Detect CRs against average and avoid objects + crcombine - Combine multiple exposures to eliminate cosmic rays + credit - Interactively edit cosmic rays using an image display + crfix - Fix cosmic rays in images using cosmic ray masks + crgrow - Grow cosmic rays in cosmic ray masks + crmedian - Detect and replace cosmic rays with median filter + crnebula - Detect and replace cosmic rays in nebular data + + overview - Overview of the package diff --git a/noao/imred/crutil/crutil.par b/noao/imred/crutil/crutil.par new file mode 100644 index 00000000..6782c988 --- /dev/null +++ b/noao/imred/crutil/crutil.par @@ -0,0 +1,3 @@ +# CRUTIL package parameter file + +version,s,h,"V1.5: August 22, 2001" diff --git a/noao/imred/crutil/doc/cosmicrays.hlp b/noao/imred/crutil/doc/cosmicrays.hlp new file mode 100644 index 00000000..55a284da --- /dev/null +++ b/noao/imred/crutil/doc/cosmicrays.hlp @@ -0,0 +1,306 @@ +.help cosmicrays Apr98 noao.imred.crutil +.ih +NAME +cosmicrays -- detect and replace cosmic rays +.ih +USAGE +cosmicrays input output +.ih +PARAMETERS +.ls input +List of input images in which to detect cosmic rays. +.le +.ls output +List of output images in which the detected cosmic rays will be replaced +by an average of neighboring pixels. If the output image name differs +from the input image name then a copy of the input image is made with +the detected cosmic rays replaced. If no output images are specified +then the input images are modified in place. In place modification of +an input image also occurs when the output image name is the same as +the input image name. +.le +.ls crmasks = "" +List of cosmic ray mask files to be created; one for each input image. If no +file names are given then no cosmic ray mask is created. If an existing +mask is specified the newly detected cosmic rays will be added. +.le + +.ls threshold = 25. +Detection threshold above the mean of the surrounding pixels for cosmic +rays. The threshold will depend on the noise characteristics of the +image and how weak the cosmic rays may be for detection. A typical value +is 5 or more times the sigma of the background. +.le +.ls fluxratio = 2. +The ratio (as a percent) of the mean neighboring pixel flux to the candidate +cosmic ray pixel for rejection. The value depends on the seeing and the +characteristics of the cosmic rays. Typical values are in the range +2 to 10 percent. This value may be reset interactively from a plot +or defined by identifying selected objects as stars or cosmic rays. +.le +.ls npasses = 5 +Number of cosmic ray detection passes. Since only the locally strongest +pixel is considered a cosmic ray, multiple detection passes are needed to +detect and replace cosmic ray events with multiple neighboring pixels. +.le +.ls window = 5 +Size of cosmic ray detection window. A square window of either 5 by 5 or +7 by 7 is used to detect cosmic rays. The smaller window allows detection +in the presence of greater background gradients but is less sensitive at +discriminating multiple event cosmic rays from stars. It is also marginally +faster. +.le +.ls interactive = yes +Examine parameters interactively? A plot of the mean flux within the +detection window (x100) vs the flux ratio (x100) is plotted and the user may +set the flux ratio threshold, delete and undelete specific events, and +examine specific events. This is useful for new data in which one is +uncertain of an appropriate flux ratio threshold. Once determined the +task need not be used interactively. +.le +.ls train = no +Define the flux ratio threshold by using a set of objects identified +as stars (or other astronomical objects) or cosmic rays? +.le +.ls objects = "" +Cursor list of coordinates of training objects. If null (the null string "") +then the image display cursor will be read. The user is responsible for first +displaying the image. Otherwise a file containing cursor coordinates +may be given. The format of the cursor file is "x y wcs key" where +x and y are the pixel coordinates, wcs is an arbitrary number such as 1, +and key may be 's' for star or 'c' for cosmic ray. +.le +.ls savefile = "" +File to save (by appending) the training object coordinates. This is of +use when the objects are identified using the image display cursor. The +saved file can then be input as the object cursor list for repeating the +execution. +.le +.ls plotfile +If a plot file is specified then the graph of the flux ratio (x100) vs +the mean flux (x100) is recorded as metacode. This may be spooled or examined +later. +.le +.ls graphics = "stdgraph" +Interactive graphic output device for interactive examination of the +detection parameters. +.le +.ls cursor = "" +Interactive graphics cursor input. If null the graphics display cursor +is used, otherwise a file containing cursor input may be specified. +.le +.ls answer +This parameter is used for interactive queries when processing a list of +images. The responses may be "no", "yes", "NO", or "YES". The upper case +responses permanently enable or disable the interactive review while +the lower case reponses allow selective examination of certain input +images. \fIThis parameter should not be specified on the command line. +If it is then the value will be ignored and the task will act as if +the answer "yes" is given for each image; i.e. it will enter the interactive +phase without prompting.\fR +.le +.ih +IMAGE CURSOR COMMANDS + +.nf +? 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 +.fi + +GRAPHICS CURSOR COMMANDS + +.nf +? Help +a Toggle between showing all candidates and only the training points +d Mark candidate for replacement (applys to '+' points) +e Mark candidates in a region for replacement (applys to '+' points) +q Quit and return to image cursor or 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 \fBgtools\fR) +<space> Print the pixel coordinates +.fi + +There are no colon commands except those for the windowing options (type +:\help or see \fBgtools\fR). +.ih +DESCRIPTION +Cosmic ray events in each input image are detected and replaced by the +average of the four neighbors. The replacement may be performed +directly on the input image if no output image is specified or if the +output image name is the same as the input image name. If a new image +is created it is a copy of the input image except for the replaced +pixels. +Optional output includes +a plot file showing the parameters of the +detected cosmic ray candidates and the flux ratio threshold used, a +cosmic ray mask identifying the cosmic rays found, and +a file of training objects marked with the image display cursor. The +cosmic ray mask may be used for display purposes, combined with other +masks, and with \fBcrfix\fR. + +This task may be applied to an image previously processed to detect +additional cosmic rays. + +The cosmic ray detection algorithm consists of the following steps. +First a pixel must be the brightest pixel within the specified +detection window (either 5x5 or 7x7). The mean flux in the surrounding +pixels with the second brightest pixel excluded (which may also be a +cosmic ray event) is computed and the candidate pixel must exceed this +mean by the amount specified by the parameter \fIthreshold\fR. A plane +is fit to the border pixels of the window and the fitted background is +subtracted. The mean flux (now background subtracted) and the ratio of +this mean to the cosmic ray candidate (the brightest pixel) are +computed. The mean flux (x100) and the ratio (x100) are recorded for +interactive examination if desired. + +Once the list of cosmic ray candidates has been created and a threshold for +the flux ratio established (by the parameter \fIfluxratio\fR, by the +"training" method, or by using the graphics cursor in the interactive plot) +the pixels with ratios below the threshold are replaced in the output by +the average of the four neighboring pixels (with the second strongest pixel +in the detection window excluded if it is one of these pixels). Additonal +pixels may then be detected and replaced in further passes as specified by +the parameter \fInpasses\fR. Note that only pixels in the vicinity of +replaced pixels need be considered in further passes. + +The division between the peaks of real objects and cosmic rays is made +based on the flux ratio between the mean flux (excluding the center +pixel and the second strongest pixel) and the candidate pixel. This +threshold depends on the point spread function and the distribution of +multiple cosmic ray events and any additional neighboring light caused +by the events. This threshold is not strongly coupled to small changes +in the data so that once it is set for a new type of image data it may +be used for similar images. To set it initially one may examine the +scatter plot of the flux ratio as a function of the mean flux. This +may be done interactively or from the optional plot file produced. + +After the initial list of cosmic ray candidates has been created and before +the final replacing cosmic rays there are two optional steps to allow +examining the candidates and setting the flux ratio threshold dividing +cosmic rays from real objects. The first optional step is define the flux +ratio boundary by reference to user specified classifications; that is +"training". To do this step the \fItrain\fR parameter must be set to yes. +The user classified objects are specified by a cursor input list. This +list can be an actual file or the image display cursor as defined by the +\fIobjects\fR parameter. The \fIsavefile\fR parameter is also used during +the training to record the objects specified. The parameter specifies a +file to append the objects selected. This is useful when the objects are +defined by interactive image cursor and does not make much sense when using +an input list. + +If the \fIobjects\fR parameter is specified as a null string then +the image display cursor will be repeatedly read until a 'q' is +entered. The user first displays the image and then when the task +reads the display cursor the cursor shape will change. The user +points at objects and types 's' for a star (or other astronomical +object) and 'c' for a cosmic ray. Note that this input is used +to search for the matching object in the cosmic ray candidate list +and so it is possible the selected object is not in the list though +it is unlikely. The selection will be quietly ignored in that case. +To exit the interactive selection of training objects type 'q'. + +If 'g' is typed a graph of all the candidates is drawn showing +"flux" vs. "flux ratio" (see below for more). Training objects will +be shown with a box and the currently set flux ratio threshold will +also be shown. Exiting the plot will return to entering more training +objects. The plot will remain and additional objects will immediately +be shown with a new box. Thus, if one wants to see the training +objects identified in the plot as one selects them from the image +display first type a 'g' to draw the initial plot. Also by switching +to the plot with 'g' allows you to draw surface plots (with 's') or +get the pixel coordinates of a candidate (the space key) to be +found in the display using the coordinate readout of the display. +Note that the display interaction is simpler than might be desired +because this task does not directly connect to the display. + +The most likely use for training is with the interactive image display. +However one may prepare an input list by other means, one example +is with \fBrimcursor\fR, and then specify the file name. The savefile +may also be used a cursor input to repeat the cosmic ray operation +(but be careful not to have the cursor input and save file be the +same file!). + +The flux ratio threshold is determined from the training objects by +finding the point with the minimum number of misclassifications +(stars as cosmic rays or cosmic rays as stars). The threshold is +set at the lowest value so that it will always go through one of +the cosmic ray objects. There should be at least one of each type +of object defined for this to work. The following option of +examining the cosmic ray candidates and parameters may still be +used to modify the derived flux ratio threshold. One last point +about the training objects is that even if some of the points +lie on the wrong side of the threshold they will remain classified +as cosmic ray or non-cosmic ray. In other words, any object +classified by the user will remain in that classification regardless +of the final flux ratio threshold. + +After the training step the user will be queried to examine the candidates +in the flux vs flux ratio plane if the \fIinteractive\fR flag is set. +Responses may be made for specific images or for all images by using +lower or upper case answers respectively. When the parameters are +examined interactively the user may change the flux ratio threshold +('t' key). Changes made are stored in the parameter file and, thus, +learned for further images. Pixels to be deleted are marked by crosses +and pixels which are peaks of objects are marked by pluses. The user +may explicitly delete or undelete any point if desired but this is only +for special cases near the threshold. In the future keys for +interactive display of the specific detections will be added. +Currently a surface plot of any candidate may be displayed graphically +in four 90 degree rotated views using the 's' key. Note that the +initial graph does not show all the points some of which are clearly +cosmic rays because they have negative mean flux or flux ratio. To +view all data one must rewindow the graph with the 'w' key or ":/" +commands (see \fBgtools\fR). +.ih +EXAMPLES +1. To replace cosmic rays in a set of images ccd* without training: + +.nf + cl> cosmicrays ccd* new//ccd* + ccd001: Examine parameters interactively? (yes): + [A scatter plot graph is made. One can adjust the threshold.] + [Looking at a few points using the 's' key can be instructive.] + [When done type 'q'.] + ccd002: Examine parameters interactively? (yes): NO + [No further interactive examination is done.] +.fi + +After cleaning one typically displays the images and possibly blinks them. +A difference image or mask image may also be created. + +2. To use the interactive training method for setting the flux ratio threshold: + +.nf + # First display the image. + cl> display ccd001 1 + z1 = 123.45 z2= 543.21 + cl> cosmicrays ccd001 ccd001cr train+ + [After the cosmic ray candidates are found the image display + [cursor will be activated. Mark a cosmic ray with 'c' and + [a star with 's'. Type 'g' to get a plot showing the two + [points with boxes. Type 'q' to go back to the image display. + [As each new object is marked a box will appear in the plot and + [the threshold may change. To find the location of an object + [seen in the plot use 'g' to go to the graph, space key to find + [the pixel coordinates, 'q' to go back to the image display, + [and the image display coordinate box to find the object. + [When done with the training type 'q'. + ccd001: Examine parameters interactively? (yes): no +.fi + +3. To create a mask image a bad pixel file must be specified. + +.nf + cl> cosmicrays ccd001 ccd001 crmask=crccd001 +.fi +.ih +SEE ALSO +crmedian, crnebula, crgrow, crfix, credit, gtools, imedit, rimcursor +.endhelp diff --git a/noao/imred/crutil/doc/craverage.hlp b/noao/imred/crutil/doc/craverage.hlp new file mode 100644 index 00000000..bd4ef7c9 --- /dev/null +++ b/noao/imred/crutil/doc/craverage.hlp @@ -0,0 +1,232 @@ +.help craverage Apr98 noao.imred.crutil +.ih +NAME +craverage -- detect CRs and objects using average filter +.ih +SYNOPSIS +\fBCraverage\fR detects cosmic rays and objects using a moving block +average filter with the central pixel plus some number of additional high +pixels excluded and a median of an annulus around the block average box. +It avoids identification of the cores of objects as cosmic rays by +excluding pixels within the detected objects as cosmic ray candidates. +.ih +USAGE +.nf +craverage input output +.fi +.ih +PARAMETERS +.ls input +List of input images in which to detect cosmic rays and objects. +.le +.ls output +List of output images in which cosmic rays are replaced by the block average +value excluding the cosmic ray. If no output image name is given then +no output image will be created. +.le +.ls crmask = "" +List of input and output cosmic ray and object masks. If the mask exists +then the mask values are used to exclude data pixels from the calculations +and zero mask values are candidates for cosmic rays or objects. +Detected cosmic rays and objects are identified in the mask with values +given by the \fIcrval\fR and \fIobjval\fR parameters. If no output cosmic +ray mask is given then no mask will be created. +.le +.ls average = "" +List of output block average filtered images. If no image name is given +then no image will be created. +.le +.ls sigma = "" +List of output sigma images. If no image name is given then no image +will be created. +.le + +.ls navg = 5 (minimum of 3) +Square block average filter size given as the number of pixels along an +edge. The value will be rounded up to an odd value to be symmetrical +around the center pixel excluded from the average. +.le +.ls nrej = 0 (minimum of 0) +Number of additional highest pixels to exclude, in addition to the +central pixel, in the block average. The value should be small but it +is needed to deal with cosmic rays that are bigger than a single pixel. +.le +.ls nbkg = 5 (minimum of 1) +Background annulus width around the box average filter in pixels. The +median of the pixels in this annulus is used to estimate the background. +.le +.ls nsig = 25 (minimum of 10) +Square box size for empirical sigma estimates given as the number of +pixels along an edge. The sigma is estimated using percentile points +of the pixels in the box. The size of the box should contain +of order 100 pixels or more. +.le +.ls var0 = 0., var1 = 0., var2 = 0. +Variance coefficients for the variance model. The variance model is + +.nf + variance = var0 + var1 * data + var2 * data^2 +.fi + +where data is the maximum of zero and the average filtered pixel value and +the variance is in data numbers. All the coefficients must be positive or +zero. If they are all zero then empirical data sigmas are estimated by a +percentile method in boxes of size given by \fInsig\fR. +.le + +.ls crval = 1 +Mask value for detected cosmic rays. It is legal for the value to be +zero to not mark the cosmic rays in the output mask. +.le +.ls lcrsig = 10., hcrsig = 5. +Low and high sigma factors for detecting cosmic rays. These factors +multiply the computed or estimated sigma at each pixel and these threshold +values are compared to the difference between the candidate pixel and the +block average filter value (average of box around the pixel). This only +applies to pixels where the block average filter value is within a +specified threshold of the background estimate; i.e. the average value is +not considered as part of an object. +.le +.ls crgrow = 0. +Cosmic ray growing radius. Pixels detected and marked in the output cosmic +ray mask by the \fIcrval\fR value are increased in size in the mask (but +not replaced in the output image) by also flagging all zero valued mask +pixels within this specified radius with the cosmic ray mask value. This +is done after the detection phase is complete. The separation between +pixels is the distance between pixel centers computed as a real value. +Note a value of at least one is required to affect other mask pixels. +.le + +.ls objval = 0 +Mask value for detected objects. It is legal for the value to be +zero to not mark the objects in the output mask. +.le +.ls lobjsig = 10., hobjsig = 5. +Low and high sigma factors for detecting objects. These factors multiply +the computed or estimated sigma at each pixel and these threshold values +are compared to the difference between the block average filter value and +the background annulus median. If the values are made very large then +object detection can be eliminated and cosmic rays will be detected +everywhere. +.le +.ls objgrow = 0. +Object detection growing radius. Pixels detected and marked in the output +mask by the \fIobjval\fR value are increased in size in the mask by also +flagging all zero valued mask pixels within this specified radius with the +cosmic ray mask value. This is done after the detection phase is complete +and so object grown pixels are not used in excluding cosmic ray +candidates. The separation between pixels is the distance between pixel +centers computed as a real value. Note a value of at least one is +required to affect other mask pixels. +.le +.ih +DESCRIPTION +\fBCraverage\fR detects cosmic rays and objects using a moving block +average filter with the central pixel and a specified number of additional +highest pixels excluded and a median of an annulus around the block average +box. It avoids identification of the cores of objects as cosmic rays by +excluding pixels within the detected objects as cosmic ray candidates. + +The block average filter computes the average of pixels in a box with the +central or target pixel excluded. In addition the \fInrej\fR parameter can +be used to exclude that number of highest remaining pixels as possible +contamination from cosmic rays which are larger than one pixel or possibly +a very nearby additional cosmic ray. The \fInrej\fR value should be kept +small relative to the total number of pixels in the average so that the +average will still be elevated over the median in real underlying objects. +The resulting average is used as the prediction for the value of the target +pixel. The median of the pixels in a square background annulus around the +block average box provides the prediction for the background at the target +pixel. + +The target pixel is considered part of an object if the difference between +the average value and the median background exceeds a specified threshold. +If the pixel is NOT considered to be part of an object then if the +difference between the pixel value and the average value exceeds a +different specified threshold it is identified as a cosmic ray. + +The thresholds are defined in terms of sigma factors, which may be +different for positive and negative deviations and for object and +cosmic ray identification. The sigma factors multiply an estimate +for the statistical sigma of the target pixel. The estimate is +either based on a noise model or sigma of pixels in a box near the +target pixel. + +The \fIcrmask\fR parameter specifies a pixel mask for the image. If the +mask exists then non-zero mask values will be used to exclude pixels from +the average, background median, and empirical sigma estimates. Also any +pixels with non-zero mask values will not be altered either in the output +image or in the final mask. If the mask does not exist then it behaves as +if all mask values are zero. If all pixels in the average box or median +annulus are previously flagged then the estimates will be undefined and +nothing will be done to the output image or mask. Because the task can +use an input mask to mark pixels not to be considered it can be used +in an iterative fashion. + +The noise model is given by the formula + +.nf + variance = var0 + var1 * data + var2 * data^2 +.fi + +where data is the maximum of zero and the average estimate for the target +pixel. The coefficients are all given in terms of the data numbers. This +model can be related to common detector parameters. For CCDs var0 is the +readout noise expressed as a variance in data numbers and var1 is the +inverse gain (DN/electrons). The second order coefficient has the +interpretation of flat field introduced variance. + +If all the coefficients are zero then an empirical sigma is estimated as +follows. The input image is divided into square blocks of size +\fInsig\fR. The (unmasked) pixel values in a block are sorted and the +pixel values nearest the 15.9 and 84.1 percentiles are selected. These are +the one sigma points in a Gaussian distribution. The sigma estimate is the +difference of these two values divided by two. This algorithm is used to +avoid contamination of the sigma estimate by the bad pixel values. The +block size must be at least 10 pixels in each dimension to provide +sufficient pixels for a good estimate of the percentile points. The sigma +estimate for a pixel is the sigma from the nearest block. A moving box is +not used for reasons of efficiency. + +If an output image name is specified then the output image is produced as a +copy of the input image but with the identified cosmic ray pixels replaced +by the average predicted value. Other optional output images are +the average filtered values and the sigma values. + +If a mask is specified the detected cosmic rays will be identified with +values given by the \fIcrval\fR parameter and object pixels will be +identified with values given by the \fIobjval\fR parameter. Note that one +does not need to use an output image and the cosmic rays can be replaced by +interpolation in the data using the tasks \fIcrfix\fR, \fIfixpix\fR, or +\fIccdproc\fR. + +One final step may be applied to the output mask. The mask values +identified with the \fIcrval\fR and \fIobjval\fR values may be grown +by identifying pixel values within a specified radius with the same +mask value. Note that this step is done at the end and so any pixels +in a preexisting input mask with the same values will also be grown. +Also the grown pixels will not affect the output cosmic ray replaced +image. See \fIcrgrow\fR for a further discussion. +.ih +EXAMPLES +This example illustrates using the \fBcraverage\fR task to +create a mask with cosmic rays and objects identified and displayed. +The image is a CCD image with a readout noise of 5 electrons +and a gain of 3 electrons per data number. This implies variance +model coefficients of + +.nf + var0 = (5/3)^2 = 2.78 + var1 = 1/3 = 0.34 +.fi + +.nf + cl> display obj001 1 # Display in first frame + cl> craverage obj001 "" crmask=mask001 var0=2.78 var1=0.34\ + >>> crval=1 objval=2 + cl> display crobj001 2 overlay=mask001 ocol="1=green,2=red" +.fi +.ih +SEE ALSO +cosmicrays, crnebula, median, crfix, crgrow, crmedian +.endhelp diff --git a/noao/imred/crutil/doc/crcombine.hlp b/noao/imred/crutil/doc/crcombine.hlp new file mode 100644 index 00000000..3dddaebe --- /dev/null +++ b/noao/imred/crutil/doc/crcombine.hlp @@ -0,0 +1,35 @@ +.help crcombine Apr98 noao.imred.crutil +.ih +NAME +crcombine -- combine multiple exposures to eliminate cosmic rays +.ih +USAGE +.nf +crcombine input output +.fi +.ih +PARAMETERS +See parameters for \fBimcombine\fR. +.ih +DESCRIPTION +This task is a version of \fBimcombine\fR. See the help for that task +for a description of the parameters and algorithms. + +For the purpose of removing cosmic rays the most useful options +are the "crreject" algorithm and/or combining with a median. Many other +options may work as well. The best use of this task depends on the +number of images available. If there are more than a few images the +images should be combined with an "average" and using a rejection +algorithm. +.ih +EXAMPLES +1. To combine two images using the gain and read noise parameters in +the image header: + +.nf + cl> crcombine obj012,obj013 abc gain=gain rdnoise=rdnoise +.fi +.ih +SEE ALSO +imcombine +.endhelp diff --git a/noao/imred/crutil/doc/credit.hlp b/noao/imred/crutil/doc/credit.hlp new file mode 100644 index 00000000..f038e5c1 --- /dev/null +++ b/noao/imred/crutil/doc/credit.hlp @@ -0,0 +1,39 @@ +.help credit Apr98 noao.imred.crutil +.ih +NAME +credit -- interactively edit cosmic rays using an image display +.ih +USAGE +.nf +credit input output +.fi +.ih +PARAMETERS +See parameters for \fBimedit\fR. +.ih +DESCRIPTION +This task is a version of \fBimedit\fR. See the help for that task +for a description of the parameters and algorithms. + +For the purpose of editing cosmic rays the most useful editing option +is 'b' to replace cosmic rays in a circular annulus using local sky +values. This can be done interactively or using a list of positions +along with the \fIdefault\fR parameter value. +.ih +EXAMPLES +1. To replace cosmic rays interactively. + +.nf + cl> credit obj012 crobj012 crmask012 +.fi + +2. To use a two column list of positions and remove the cosmic rays using +the 'b' key algorithm. + +.nf + cl> credit obj012 crobj012 cursor=crlist.dat display- +.fi +.ih +SEE ALSO +imedit, epix +.endhelp diff --git a/noao/imred/crutil/doc/crfix.hlp b/noao/imred/crutil/doc/crfix.hlp new file mode 100644 index 00000000..5794a001 --- /dev/null +++ b/noao/imred/crutil/doc/crfix.hlp @@ -0,0 +1,48 @@ +.help crfix Apr98 noao.imred.crutil +.ih +NAME +crfix -- fix cosmic rays in images using cosmic ray masks +.ih +USAGE +.nf +crfix input output masks +.fi +.ih +PARAMETERS +.ls input +Input two dimensional image to be "fixed" (modified) by linear interpolation. +.le +.ls output +Output image. If the output image name exactly matches the input +image name (including extensions) then the image will be modified in place. +.le +.ls crmask +Cosmic ray mask identifying the cosmic rays to be fixed. The mask +values are zero for good data and non-zero for cosmic rays. +.le +.ih +DESCRIPTION +The input and output images are specified by the \fIinput\fR and +\fIoutput\fR parameters. If the input and output image names are +identifical (including extensions) then image is modified in place. Cosmic +rays, identified in a cosmic ray mask specified by the \fIcrmask\fR +parameter, are replaced in the output image by linear interpolation along +lines or columns using the nearest good pixels. The special mask name +"BPM" may be used to select a mask name given in the input image header +under the keyword "BPM". + +Cosmic ray pixels are "fixed" by replacing them with values +by linear interpolation to the nearest pixels not identified as bad. +The interpolation direction is the shortest length between good pixels +along columns or lines. +.ih +EXAMPLES +1. To replace cosmic rays in an image: + +.nf + cl> crfix obj012 crobj012 crmask012 +.fi +.ih +SEE ALSO +fixpix, crmedian, crnebula, cosmicrays, credit, epix +.endhelp diff --git a/noao/imred/crutil/doc/crgrow.hlp b/noao/imred/crutil/doc/crgrow.hlp new file mode 100644 index 00000000..f3a1ff5c --- /dev/null +++ b/noao/imred/crutil/doc/crgrow.hlp @@ -0,0 +1,55 @@ +.help crgrow Apr98 noao.imred.crutil +.ih +NAME +crgrow -- grow cosmic rays in cosmic ray masks +.ih +USAGE +crgrow input output radius +.ih +PARAMETERS +.ls input +List of cosmic ray masks to be modified. +.le +.ls output +List of output modified cosmic ray masks. The input and output lists must +match. If the input and output cosmic ray masks are specified as the same +then the input mask will be modified in place. +.le +.ls radius = 1. +Replacement radius around cosmic rays. +If a pixel is within this distance of a cosmic ray pixel +it is identified by a value of 1 in the output cosmic ray mask. Distances are +measured between pixel centers which are have integer coordinates. +.le +.ls inval = INDEF +Mask value to be grown. A value of INDEF will grow all non-zero values. +.le +.ls outval = INDEF +Mask value for grown pixels. A value of INDEF will use the value of the +pixel being grown for the grown pixel value. +.le +.ih +DESCRIPTION +The cosmic ray pixels, identified by the "inval" parameter, in the input +mask are located and all unmasked (zero valued) pixels within the specified +grow radius are set to a value given by the "outval" parameter. The +distance between pixels is measured as a cartisian logical pixel coordinate +distance. +.ih +EXAMPLES +1. A radius of 1 will grow cosmic rays in a "plus" pattern. + +.nf + cl> crgrow crmask1 crmask2 1 +.fi + +2. A radius of 1.5 will grow cosmic rays in a box pattern. The following +will modify the input mask. + +.nf + cl> crgrow crmask crmask 1.5 +.fi +.ih +SEE ALSO +imreplace +.endhelp diff --git a/noao/imred/crutil/doc/crmedian.hlp b/noao/imred/crutil/doc/crmedian.hlp new file mode 100644 index 00000000..3147b676 --- /dev/null +++ b/noao/imred/crutil/doc/crmedian.hlp @@ -0,0 +1,157 @@ +.help crmedian Apr98 noao.imred.crutil +.ih +NAME +crmedian -- detect, fix, and flag cosmic rays using median filtering +.ih +USAGE +.nf +crmedian input output +.fi +.ih +PARAMETERS +.ls input +Input image in which to detect cosmic rays. +.le +.ls output +Output image in which cosmic rays are replaced by the median value. +If no output image name is given then no output image will be created. +.le +.ls crmask = "" +Output cosmic ray mask. Detected cosmic rays (and other deviant pixels) +are identified in the mask with values of one and good pixels with a values +of zero. If no output cosmic ray mask is given then no mask will be +created. +.le +.ls median = "" +Output median filtered image. If no image name is given then no output will be +created. +.le +.ls sigma = "" +Output sigma image. If no image name is given then no output will be +created. +.le +.ls residual = "" +Output residual image. This is the input image minus the median filtered +image divided by the sigma image. Thresholds in this image determine the +cosmic rays detected. If no image name is given then no output will be +created. +.le +.ls var0 = 0., var1 = 0., var2 = 0. +Variance coefficients for the variance model. The variance model is + +.nf + variance = var0 + var1 * data + var2 * data^2 +.fi + +where data is the maximum of zero and median pixel value and the variance +is in data numbers. All the coefficients must be positive or zero. If +they are all zero then empirical data sigmas are estimated by a percentile +method in boxes of size given by \fIncsig\fR and \fInlsig\fR. +.le +.ls lsigma = 10, hsigma = 3 +Positive sigma factors to use for selecting pixels below and above +the median level based on the local percentile sigma. Cosmic rays will +appear above the median level. +.le +.ls ncmed = 5, nlmed = 5 +The column and line size of a moving median rectangle used to estimate the +uncontaminated local image. +.le +.ls ncsig = 25, nlsig = 25 +The column and line size of regions used to estimate the uncontaminated +local sigma using a percentile. The size of the box should contain +of order 100 pixels or more. +.le +.ih +DESCRIPTION +\fBCrmedian\fR detects cosmic rays from pixels deviating by a specified +statistical amount from the median at each pixel. It outputs and set of +the following: a copy of the input image with cosmic rays replaced by the +median value, a cosmic ray mask identifying the cosmic rays, the median +filtered image, a sigma image where each pixel has the estimated sigma, and +the residual image used in detecting the cosmic rays. + +The residual image is computed by subtracting a median filtered version +of the input data from the unfiltered input data and dividing by an +estimate of the pixel sigmas. The median filter +box size is given by the \fIncmed\fR and \fInlmed\fR parameters. +If a name for the median image is specified the median filtered image +will be output. The variance at each pixel is determined either from +a variance model or empirically. If a name for the sigma image is specified +then the sigma values (the square root of the variance) will be output. +If a name for the residual image is given then the residual image +will be output. + +The empirical variance model is given by the formula + +.nf + variance = var0 + var1 * data + var2 * data^2 +.fi + +where data is the maximum of zero and median pixel value and the variance +is in data numbers. This model can be related to common detector +parameters. For CCDs var0 is the readout noise expressed as a variance in +data numbers and var1 is the inverse gain (DN/electrons). The second order +coefficient has the interpretation of flat field introduced variance. + +If all the coefficients are zero then an empirical sigma is estimated +as follows. The input image is divided into blocks of size +\fIncsig\fR and \fInlsig\fR. The pixel values in a block are sorted +and the pixel values nearest the 15.9 and 84.1 percentiles are +selected. These are the one sigma points in a Gaussian distribution. +The sigma estimate is the difference of these two values divided by +two. This algorithm is used to avoid contamination of the sigma +estimate by the bad pixel values. The block size must be at least 10 +pixels in each dimension to provide sufficient pixels for a good estimate +of the percentile points. The sigma estimate for a pixel is the sigma +from the nearest block. A moving box is not used for efficiency. + +The residual image is divided by the sigma estimate at each pixel. +Cosmic rays are identified by finding those pixels in the +residual image which have values greater than \fIhsigma\fR and bad +pixels with values below \fIlsigma\fR are also identified. + +If an output image name is specified then the output image is produced as a +copy of the input image but with the identified cosmic ray pixels replaced +by the median value. If an output cosmic ray mask is specified a cosmic +ray mask will be produced with values of zero for good pixels and one for +bad pixels. The cosmic ray mask is used to display the cosmic ray +positions found and the cosmic rays can be replaced by interpolation (as +opposed to the median value) using the task \fIcrfix\fR. + +The \fBcrmedian\fR detections are very simple and do not take into account +real structure with scales of a pixel. Thus this may clip the cores of +stars and narrow nebular features in the data. More sophisticated +algorithms are found in \fBcosmicrays\fR, \fIcraverage\fR, and +\fBcrnebula\fR. The median, sigma, and residual images are available as +output to evaluate the various aspects of the algorithm. +.ih +EXAMPLES +This example illustrates using the \fBcrmedian\fR task to +give a cosmic ray removed image and examining the results with an image +display. The image is a CCD image with a readout noise of 5 electrons +and a gain of 3 electrons per data number. This implies variance +model coefficients of + +.nf + var0 = (5/3)^2 = 2.78 + var1 = 1/3 = 0.34 +.fi + +.nf + cl> display obj001 1 # Display in first frame + cl> # Determine output image, cosmic ray mask, and residual image + cl> crmedian obj001 crobj001 crmask=mask001 resid=res001\ + >>> var0=2.78 var1=0.34 + cl> display crobj001 2 # Display final image + cl> display mask001 3 zs- zr- z1=-1 z2=2 # Display mask + cl> display res001 4 zs- zr- z1=-5 z2=5 # Display residuals +.fi + +By looking at the residual image the sigma clippig threshold can be +adjusted and the noise parameters can be tweaked to minimize clipping +of real extended structure. +.ih +SEE ALSO +cosmicrays, craverage, crnebula, median, crfix, crgrow +.endhelp diff --git a/noao/imred/crutil/doc/crnebula.hlp b/noao/imred/crutil/doc/crnebula.hlp new file mode 100644 index 00000000..174a6771 --- /dev/null +++ b/noao/imred/crutil/doc/crnebula.hlp @@ -0,0 +1,139 @@ +.help crnebula Apr98 noao.imred.crutil +.ih +NAME +crnebula -- create a cosmic ray mask from nebular images +.ih +USAGE +.nf +crnebula input output +.fi +.ih +PARAMETERS +.ls input +Input image in which cosmic rays are to be detected. +.le +.ls output +Output image in which cosmic rays are to be replaced by the median. +If no output image is given (specified as "") then no output image +is created. +.le +.ls crmask = "" +Output cosmic ray mask identifying the cosmic rays found. The mask +will have values of one for cosmic rays and zero for non-cosmic rays. +If no output cosmic ray mask is given (specified as "") then no mask +is created. +.le +.ls residual = "" +Output residual image. This is the input image minus the median filtered +image divided by the estimated sigma at each pixel. Thresholds in this +image determine the cosmic rays detected. If no image name is given then +no output will be created. +.le +.ls rmedresid = "" +Output image for the difference between the box median filter image and +the ring median filtered image divided by the estimated sigma at each +pixel. If no image name is given then no output will be created. +.le +.ls var0 = 0., var1 = 0., var2 = 0. +Variance coefficients for the variance model. The variance model is + +.nf + variance = var0 + var1 * data + var2 * data^2 +.fi + +where data is the maximum of zero and median pixel value and the variance is in +data numbers. All the coefficients must be positive or zero. If they are +all zero then empirical data sigmas are estimated by a percentile method in +boxes of size given by \fIncsig\fR and \fInlsig\fR. +.le +.ls sigmed = 3. +Sigma clipping factor for the residual image. +.le +.ls sigdiff = 3. +Sigma clipping factor for the residuals between the box median and ring median +filtered images. +.le +.ls mbox = 5 +Box size, in pixels, for the box median filtering. +.le +.ls rin = 1.5, rout = 6. +Inner and outer radii, in pixels, for the ring median filtering. +.le +.ls verbose = no +Print some progress information? +.le +.ih +DESCRIPTION +This task uses a combination of box median filtering to detect cosmic rays +and the difference between box and ring median filtering to identify +regions of fine nebular structure which should not be treated as cosmic +rays. The output consists of some set of the input image with cosmic rays +replaced by the median, a cosmic ray mask, the residual image used to +detect the cosmic rays, and the residual image used to exclude cosmic rays +in regions of nebular fine structure. The cosmic ray mask may be used +later with \fBcrgrow\fR and \fBcrfix\fR to grow and remove the cosmic rays +from the data by interpolation rather than the median. + +The algorithm is as follows. The input image is median filtered using a +box of size given by \fImbox\fR. The residual image between the unfiltered +and filter data is computed. The residuals are divided by the estimated +sigma of the pixel. Cosmic rays are those which are more than \fIsigmed\fR +above zero in the residual image. This residual image may be output if an +output name is specified. This part of the algorithm is identical to that +of the task \fIcrmedian\fR and, in fact, that task is used. + +The median image not only enhances cosmic rays it also enhances narrow fine +structure in the input image. To avoid identifying this structure as +cosmic rays a second filtered residual image is created which +preferentially identifies this structure over the cosmic rays. The input +image is filtered using a ring median of specified inner and outer radius. +The inner radius is slightly larger than the scale of the cosmic rays and +the outer radius is comparable to the box size of the box median filter. A +ring filter replaces the center of the ring by the median of the ring. The +difference between the input and ring median filtered image divided by the +estimated sigma will then be very similar to the box median residual image both +where there are cosmic rays and where there is diffuse structure but will +be different where there are linear fine structure patterns. The +difference between the median residual image and this ring median residual +image highlights the regions of fine structure. If a image name is specified +for the difference of the residual images it will be output. + +The difference of the median residual images is used to exclude any cosmic +ray candidate pixels determined from sigma clipping the box median residual +image which lie where the difference of the median residual images is +greater than \fIsigdiff\fR different from zero (both positive or +negative). + +To understand this algorithm it is recommended that the user save the +residual and residual difference images and display them and blink against +the original data. +.ih +EXAMPLES +This example, the same as in \fBcrmedian\fR, illustrates using the +\fBcrnebual\fR task to give a cosmic ray removed image and examining the +results with an image display. The image is a CCD image with a readout +noise of 5 electrons and a gain of 3 electrons per data number. This +implies variance model coefficients of + +.nf + var0 = (5/3)^2 = 2.78 + var1 = 1/3 = 0.34 +.fi + +.nf + cl> display obj001 1 # Display in first frame + cl> # Determine output image, cosmic ray mask, and residual images + cl> crnebula obj001 crobj001 crmask=mask001 resid=res001\ + >>> rmedresid=rmed001 var0=2.78 var1=0.34 + cl> display crobj001 2 # Display final image + cl> display res001 3 zs- zr- z1=-5 z2=5 # Display residuals + cl> display rmed001 4 zs- zr- z1=-5 z2=5 +.fi + +By looking at the residual image the sigma clippig threshold can be +adjusted and the noise parameters can be tweaked to minimize clipping +of real extended structure. +.ih +SEE ALSO +cosmicrays, crmedian, median, rmedian, crfix, crgrow +.endhelp diff --git a/noao/imred/crutil/doc/overview.hlp b/noao/imred/crutil/doc/overview.hlp new file mode 100644 index 00000000..cb4dc3de --- /dev/null +++ b/noao/imred/crutil/doc/overview.hlp @@ -0,0 +1,76 @@ +.help overview Apr98 noao.imred.crutil + +.ce +\fBThe Cosmic Ray Package: CRUTIL\fR + +The cosmic ray package provides tools for identifying and removing cosmic +rays in images. The tasks are: + +.nf + cosmicrays - Remove cosmic rays using flux ratio algorithm + craverage - Detect CRs against average and avoid objects + crcombine - Combine multiple exposures to eliminate cosmic rays + credit - Interactively edit cosmic rays using an image display + crfix - Fix cosmic rays in images using cosmic ray masks + crgrow - Grow cosmic rays in cosmic ray masks + crmedian - Detect and replace cosmic rays with median filter + crnebula - Detect and replace cosmic rays in nebular data +.fi + +The best way to remove cosmic rays is using multiple exposures of the same +field. When this is done the task \fBcrcombine\fR is used to combine the +exposures into a final single image with cosmic rays removed. The images +are scaled (if necessary) to a common data level either by multiplicative +scaling, an additive background offset, or some combination of both. +Cosmic rays are then found as pixels which differ by some statistical +amount away for the average or median of the data. + +A median is the simplest way to remove cosmic rays. This is an option +with \fBcrcombine\fR. But this does not make optimal use of the data. +An average of the pixels remaining after some rejection operation is better. +If the noise characteristics of the data can be described by a gain and +read noise then cosmic rays can be optimally rejected using the +"crreject" algorithm. This works on two or more images. There are +a number of other rejection algorithms which can be used as described in +the task help. + +The rest of the tasks in the package are used when only a single exposure +is available. These include interactive editing with \fBcredit\fR. The +replacement algorithms in this task may also be used non-interactively if +you have a list of pixel coordinates as input. Other tasks automatically +identifying pixels which are significantly higher than surrounding pixels. + +The simplest of these tasks is \fBcrmedian\fR. This replaces +cosmic rays with a median value and produces a cosmic ray +mask which is a simple type of integer image where good pixels have a value +of zero and bad pixels have a non-zero value. The tasks \fBcrgrow\fR and +\fBcrfix\fR are provided to use this type of cosmic ray mask. The former +will flag additional pixels within some radius of the flagged pixels in the +mask. The latter is the basic tool for replacing the identified pixels in +the data by neighboring data. It uses linear interpolation along lines or +columns. The median task is simple but it often will flag the cores of +stars or other small but real features. + +The task \fBcraverage\fR is similar to \fBcrmedian\fR in that it compares +the pixel values against a smoothed version. Instead of a median it uses +an average with the central pixel excluded. It is more sophisticated +in that it also compares the average against a larger median to see if +the region corresponds to an object. Thus it can detect objects and +the task could be used as a simple object detection task in its own right. +Because the hardest part of cosmic ray detection from a single image is +avoiding truncation of the cores of stars this task does not allow cosmic +rays to be detected where it thinks there is an object. This task is +also more versatile in allow separate mask values and works on a list +of images. + +Somewhat more sophisticated algorithms are available in the tasks +\fBcosmicrays\fR and \fBcrnebula\fR. These attempt to determine if a +deviant pixel is the core of a star or part of a linear nebular feature +respectively. + +The best use of these tasks is to experiment and iterate. In particular, +one may want to iterate a task several times and use both \fBcosmicrays\fR +and \fBcraverage\fR. + +Good hunting! +.endhelp diff --git a/noao/imred/crutil/mkpkg b/noao/imred/crutil/mkpkg new file mode 100644 index 00000000..3a55a03a --- /dev/null +++ b/noao/imred/crutil/mkpkg @@ -0,0 +1,8 @@ +# Make the package. + +$call update@src +$exit + +update: + $call update@src + ; 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 |