Program: IMAGE Purpose: imager and/or destriper for IRAS Category: IRAS File: image.src Author: Do Kester Keywords: (In the order in which they are asked for the first time) *** LOGLEVEL= amount of output in logfile [Quiet] Silent : no output at all Quiet : only global output Normal : global plus last iteration Loud : global plus all iterations IRSET= Name(s) of input IR data set(s) [quit] CENTER= Long and Lat of the map [IRDS-center] SIZE= Long and Lat size of map [IRDS-size] COOR= Coordinate system [IRDS-coor] NAME= name of the map [none] MEDIAN= apply median filter else a simple mean [no] BEAM= Radius of convolution beam in pixels [0.0] (ie. true size) PROJ= Projection type [gnomon] *** PRCEN= Long and Lat of projection center [center] *** PARSTATUS= Start and end status of destripe parameters. [new,delete] old : continue upon params from a previous session. new : start anew, irrespective of previous sessions. keep : keep the params for further use. delete: delete params at all levels. *** FPROB= Probabilty level for accepting stripe parameters [0.05] The stripe has to stand out quite well to be accepted. Use FPROB >= 1 to always accept the parameters eg. when using the destripe parameters later in HIRAS. It does no harm. SET= map to be used as template in destriping; : use map filled with values [none] ITER= number of iterations [0] PIXELS= (set of) pixel size(s) [quit] TUNE= destripe power 0: no tuning, 1: constant, 2: linear [2] ADJUST= list of snip.det's to be adjusted to map. [none] EXCLUDE= list of snip.det's to be skipped. [none] *** DETS= which detectors to use large : use large and three-quarter detectors only small : use small detectors only all : use all detectors : use the list OUTSET= (Generic name of) the output set. [none] The final map will have this name; Earlier iterations, if they are made, will be named _. OVERWRITE= Overwrite existing map. [yes] MAPNRS= list of iterations for which a map is to be made. [last] all : produce a map for every iteration. last: produce a map only for the last iteration. : produce a map for the given iterations. *** PLANES= produce plane(s) in the output map of [1] 1 : intensity map. 0 : coverage map. 2 : standard deviations map. INTENDED= Use intended positions if no bphf available [no] *** XYPLOT= give the iteration number for which to plot an overlay [0] GRDEVICE= device to plot on [?] *** FLAGS= Use also those samples flagged as [P,T] PointSource Tail of a pointsource ( only at 12 and 25 ) Glitch *** SLOW= speed of instrument wrt survey ( for AOs only ) SHOW= produce a plot with destripe results All : all snip_det will be shown without user interaction One : all snip_det will be shown after user interaction None: no display Last: display last iteration only The last mnemonic (or no mnemonic at all) can be followed by : only this combi is shown. Updates: 0.1 07 May 1991 DK. first try. 1.1 05 Dec 1991 DK: several bugs repaired 2.0 18 Mar 1992 DK, new version. 2.1 12 May 1992 DK, sloppy corners & input set option repaired 3.0 15 Jun 1992 DK. 3.1 12 Aug 1993 DK, epoch read and written as double Description: 1. Introduction keywords: IRSET= IMAGE is a GIPSY program to combine the snips of an IRDS into an image. There are several routes through the program which can be chosen to arrive at the best scientific and/or visual image. The program consists largely of two parts: one to make a map by coaddition of an IRDS and the other part uses an already made map as a template to adjust snips to, which are then coadded into a new map. This is an iterative procedure called destriping. Both parts can be run separately, although the most natural way is to coadd a map first and destripe it in the same run. In the coadd an optional median (huber) filter can be switched on. 2. Coadd keywords: CENTER= SIZE= COOR= BEAM= PIXELS= PROJ= PRCEN= NAME= A map is defined by the longitude and latitude of the center of the map, the length of the sides and the size of pixels - the pixels are always square, the map might be rectangular -, the coordinate system the coordinates are given in and the type of projection to transform spherical coordinates to a plane. Optionally a projection center different from the center of the map can be defined. To coadd the snips into a map we define two identical pixel grids. The pixels can be defined arbitrarily small. The position of each sample is projected onto the grid. Those pixels falling within a pre-chosen beam centered on the sample position are increased in the first grid with the brightness of the sample. In the other grid the corresponding pixel values are increased by one. In this way the sum of the brightnesses of the overlapping samples is stored in the first grid and the number of coverages in the second. Once all samples are processed we divide the first grid by the second, pixel by pixel, yielding a map in which every pixel represents the average of all overlapping samples. The beam can be chosen arbitrarily, provided that it is larger than the detector size and provided that at least one pixel falls within the beam for each beam position. The default value beam=0.0 chooses a convolution beam which has the same shape and size as the detectors themselves; see IRAS ExplSupp p. II-12 & 13. This yields the highest resolution possible in coaddition, although it is different in inscan and cross scan directions. Dead detectors can not be processed by IMAGE; small detectors are not processed by default. Section 5 explains how to exclude other detectors or to include the small ones. Table 1 defines dead and small detectors. =========================================================================== Table 1. Dead and small detectors band 1 band 2 band 3 band 4 dead 17 20 36 small 26 47 39 46 11 31 55 62 =========================================================================== The keywords CENTER, SIZE, PIXELS and PRCEN are angular reals. An angular real is defined as a set of reals separated by 'd', 'm', 's' or 'h'. `d' indicates that the preceeding value is in degrees. `m' indicates that the preceeding value is in minutes. `s' indicates that the preceeding value is in seconds. `h' indicates that the preceeding value is in hours; following 'm' and/or 's' reals are converted to hourly values. The sign of the first real is given to the whole value. The default letter is 'd'. Upper and lower case letters are equivalent. Examples -50d30m translates to -50.5 degrees 12h12.5s 180.0521 33.33 33.33 3. Destripe keywords: ITER= PIXELS= TUNE= PARSTATUS= The destriping procedure boils down to three steps which are taken for each detector snip. First the trajectory of the detector over the map is determined; second the difference between the map values on the trajectory and the corresponding snip values is fitted with the function the user requested and third the detector snip values, corrected with the function just determined, are coadded into a new map. When all the snips are processed in this manner, the new map resembles the old map in a global way. They are made from the same set of snips and the old map functioned as a template for the new map. Locally, however the differences can be appreciable. If the snips which constitute the old map are not exactly on the same scale and offset, the map will show stripes. The stripes being the average of all locally overlapping snips are not as large as the differences between the snips themselves. By fitting the detector snips to a (striped) map the differences between the detectors are lessened and the next map will show less stripes. Iterating the whole procedure several times yields a map in which the snip pattern is (almost) invisible. The dominant snip direction can sometimes still be perceived as structured background noise. The fit of the differences between the map values and the detector values is performed using a robust estimator (see D. Kester IPAC IOM 701-039-89(1) and references therein). The particular estimator was a Tukey, which smoothly diminishes the influence of points as the stray away from the fit and disregards them completely if they are further than 7 times a globally determined scale factor. One way to reduce the stripes in a map is to increase the pixel size. More detectors from different snips fall on the same pixel and the differences are averaged out. Although such a map is clearly not desirable as an end product due to its coarse grain, it can be used as template to fit the snips to, which are then coadded into a new, finer grained map. Two precautions should be taken in this scheme. The first is that one large point source might `spoil' a pixel much larger in size than the point source. Neighbouring snips which have not seen the point source are compared with the spoiled pixel, possibly resulting in an improper fit. To remedy this, point sources are detected with a very simple algorithm and excluded from the coadd of those maps that are not yet on their final grid; see section 6. Another way around problems like these is to introduce a mask; see section 7. The second danger is that a real slope present in the snips that should eventually show up in the map might not be conserved, especially in the corners of the map where a detector snip might overlap only one pixel. This is countered by limiting the fit to an offset, the most important effect of the stripes, in the coarser maps. The number of iterations, the power of the fit and the pixel size per iteration can be chosen by the user. The destripe parameters are (temporarily) stored as descriptor item in the header of the IRDS. For each detector snip a set of 3 parameters can be generated: offset, drift and standard deviation which adds another 20% to the already large descriptor file of the IRDS. The (hidden) keyword PARSTATUS determines the status of the parameters at the start and at the end of the program. It can take two values. One from the set [old,new], which selects whether one wants to continue with previously generated (old) parameters, or whether one want a new session. And one from the set [keep,delete] which determines whether the parameters should be kept for further processing or not. The default is new,delete. 4. Median Filter keywords: MEDIAN= The median filter is build into the program to eliminate transient phenomena as glitches and glints, space junk and asteroids from the map. All these phenomena are only present in a minority of the snips, but occasionally they might be (very) bright. A straight arithmetic mean leaves too much of it in the map. A median or a median-like filter should perform better. The median might be the ideal choice if not it requires that all the samples on all the pixels are stored individually until all snips are processed. Then a median has to be made by ordering the values. It requires an awfull lot of memory and a non-trivial sorting. The IMAGE program has instead a median-like (Huber) filter. It is a combination of a straight average for the central points and a median for the outliers. As explained in Kester (op.cit.) there is a simple way to calculate the huber average and scale factor which entails no storage of large numbers of values. It is an iterative procedure in its own right and it requires a standard deviations map to be made. Above it was mentioned that the transient phenomena are a minority of the sample values falling on a pixel. In a coarse grained map point sources themselves might become a minority of all the overlapping pixel values. Performing a median filter in such a situation would eliminate the point sources too. So the median filter is only applied when the map is on its final pixel grid and the user is responsible to assure the the point sources do not drown in the final map. The program does not make a standard deviations map as long as the median filter is off. The median filter can be switched on/off by the user. 5. Selection keywords: EXCLUDE= ADJUST= DETS= If for some reason one or more snips and/or detectors are deemed unreliable, improperly corrected, containing left-over zodiacal emission or whatever, they can be handled as a group in a different way. The choices are to exclude them completely from the whole procedure or to adjust them to the template made up by the other snips. Snips in the exclude list are just disregarded as if they never existed. Although a record of them can appear in the administration file (section 7). The snips in the adjust list are excluded in the first coadd of the map and in subsequent iterations they are taken along with all the others. Both EXCLUDE and ADJUST are lists of snip.det i.e. real numbers where the integral part equals the snip number and the fractional part equals the detector number divided by 100. When the integral part or the fractional part equals 0 it means that it is valid for all snips resp. all detectors. Examples: list item snip detector band 12.39 12 39 2 12.01 12 1 4 12 12 all all 0.1 all 10 3 Duplicate items and items referring to a detector which is not present in the IRDS are silently removed. The keyword DETS determines which detectors are coadded. It can take the mnemonic `large' which means that only the large and the three-quarter size detectors are used. This is the default and it is equivalent to NOT using the small detectors of table 1. The mnemonic `small' means that only the small detectors are used and `all' means that all detectors are used. In stead of mnemonics it can take a list of detectors to be used. Again, dead or duplicate detectors and detectors which are not present in the IRDS are silently removed. 6. Sample selection keyword : FLAGS= During the processing of a snip two simple filters are passed over the data stream: one to detect glitches and the other to detect point sources. The glitch filter extends over 3 points [-1,+2,-1]. The pointsource filter extends over 9 points [-1,-1,-1,+2,+2,+2,-1,-1,-1]. When the output of the point source filter exceeds a certain level, that sample and some before and after are flagged as source. How many depends on the strength of the source. If the band is 12 or 25 micron a number of points following the source are flagged as tail; again how many depends on the strength of the source. A sample is flagged as a glitch if it is not a source and the glitch filter exceeds a certain level. Samples flagged as sources, tails or glitches are not coadded into map which are not yet on their final pixel scale (see section 3). Which of these samples are coadded into the final map is under control of the user with the (hidden) keyword FLAGS. It can have 3 values at most and it checks the first letter of the mnemonics: Point Source indicating both the source flag Tail indicating the tail flag Glitch indicating the glitch flag. The values are `added'. By default the value of FLAGS=p,t indicating that in the final map both samples with point source flags and with tail flags are present. 7. Input map and mask keywords: SET= MASK= In stead of using a coadded map as template for the destriping, the user can give an already made map, maybe from another band or from another instrument, to be used as template. This input map should have the same center and size as the map to be made. It should be a GDS set. Another option is to give just a number in answer to the keyword SET=. The program will then generate a template map with pixel values equal to the number, the most obvious choice being 0. A mask is a map with some values set to (GIPSY)blank. All non-blank values are considered the same: only 1 of bit information is used. It is stored parallel to the template map and read at the same time during the destriping process. Whenever a value in the mask equals blank, the corresponding difference between snip and template is not taken along in the fit. The mask is propagated through the destripe iteration, appropriately adjusted for a change in pixel size. The GIPSY task `blot' can be used to produce a mask. 8. Output map keywords: OUTSET= OVERWRITE= MAPNRS= PLANES= The map(s) are written as GIPSY sets to a file of which the basic name is supplied by the user. The last iteration has this basic name. For earlier iterations, if they must be written, the basic name is appended with an underscore followed by the iteration number of which the map is the result. The iterations for which a map will be written are under control of the user although the last iteration is always written although once a basic name is supplied. In the 3rd dimension of the set the standard deviations map and/or the coverage map might be stacked. The former is a map which gives pixel by pixel the standard deviation of all samples which overlay a pixel in the intensity map. The latter is a map which sums the weights of all samples which fall on the pixels in the intensity map. The number and the order in the 3rd dimension of the GIPSY set can also be controlled. 9. Logfile keywords: LOGLEVEL= The administration file logs the setting of the parameters and the proceeding of the program. The loglevel `silent' provides no information about the run of the program; the level `quiet' echoes global information about the IRDS and about the overall fit of the stripes; the level `normal' provides the global information plus detailed information about the fit in the last iteration and finally, the level `loud' yields the global information plus detailed information about the fit in each iteration. Setting the GIPSY mode to `test' will give the loud information plus an unspecified amount of other garbage. 10. Graphics keywords: XYPLOT= GRDEVICE= SHOW= During the program 2 kinds of graphic output can be produced. The first shows where each of the samples fall within the map. With XYPLOT it can be chosen for which iteration it will be displayed, although the information is exactly the same in all iterations. The second kind of graphic output displays the values of the snip and the values under its cross section over the template map against time. This is in the upper panel. The difference between the two is shown in the lower panel together with the fit, derived by the program. In the middle panel it shows the timing of several flags (if any is present). Where there appears a flag, the difference data are not taken along in the fit. It is under control of the user how snips are displayed with SHOW. The answer `all' will display all snips as they are processed (auto); the answer `one' will wait after each display until a new answer to SHOW has been given (manual); `none' will give no display at all and `last' will display only the last iteration. If `last' is given or if no answer is given, a list of 3 integers can follow which are interpreted as sop, att and detector. Only that combination is displayed. One or more of these numbers can be zero which is interpreted as a wildcard. Example: answer to SHOW displays iteration sop att det auto/man L 350 23 12 last 350 23 12 manual 120 12 all 120 12 all manual 0 0 59 all all all 59 auto A all all all all auto L last all all all auto N none none none none auto O all all all all manual As it takes quite some time to set up the display, the fastest option is `none'. At any moment during the execution SHOW can be reset. NB. It is left to the user to make sure that the two kinds of graphic display are not used at the same time. 11. Memory usage. To estimate the memory necessary for this program, ie. for the internal storage of the map and the template, which are the main culprits, calculate the number of pixels in your map, multiply by 3 (if median) or 2 (if not) and multiply by 4 (bytes/real). The total is the number of bytes the program needs in two continuous arrays, one for the map and the other for the template. For a 1000 square pixels map this accumulates to about 25 Mbyte. On top of that some memory is needed e.g. to accustom one snip at a time.