;----------------------------------------------------------------------------- pro trace_prep, input1, input2, index_out, image_out, $ normalize=normalize, $ darkdir=darkdir, user_dark=user_dark, $ no_darksub=no_darksub, $ flatdir=flatdir, no_flatfield=no_flatfield, $ unspike=unspike, $ destreak=destreak, deripple=deripple, $ wave2point_correct=wave2point_correct, $ float=float, new_avg=new_avg, $ subimgx=subimgx, subimgy=subimgy, $ sllex=sllex, slley=slley, $ nodata=nodata, no_calib=no_calib, $ original=original, outminsiz=outminsiz, $ outdir=outdir, outtrfits=outtrfits, $ outflatfits=outflatfits, qstop=qstop, $ quiet=quiet, verbose=verbose, $ n_pixel=n_pixel, run_time=run_time ; *********************************************************** ; despike=despike (TBD) ; response_norm=response_norm (TBD) ;+ ; NAME: ; TRACE_PREP ; ; PURPOSE: ; Process TRACE image(s). ; ; A preliminary version of a processing routine for preparing the ; TRACE image data for analysis. The steps performed are: ; 1. Read in raw FITS image(s) from a filelist and ; decompress the data and structure; ; or read in a datacube and structure ; 2. Fill pixels of value = 0 with mean pixel value of entire ; image ; 3. Replace near saturated pixels with values > 4095 ; 4. Subtract the dark pedestal (ADC offset) & current from each ; image ; 5. Perform ccd gain calibration and lumogen degradation ; correction by dividing by time and wavelength dependent ; flat field ; 6. Options to remove radiation-belt/cosmic-ray spikes and ; streaks ; 7. Option to remove background diffraction pattern (ripple) ; 8. Option to normalize each image for exposure ; 9. Option to extract a subimage from each image ; 10. Output the corrected image(s) in an updated structure ; and data cube, and optionally output a FITS file ; with 1 or more images plus a binary extension or as ; 2D flat FITS files (1 per image) ; In the future this routine may also: ; 11. Instrument response normalization to physical units ; ; CATEGORY: ; FITS processing ; ; SAMPLE CALLING SEQUENCE: ; TRACE_PREP, input1, input2, [index_out, image_out] [,/outtrfits] ; [,outdir=outdir] [,/normalize] [,darkdir=darkdir] ; [,/user_dark] [,user_dark=user_dark] [,/no_darksub] ; [,flatdir=flatdir] [,/no_flatfield] [,unspike=unspike] ; [,/destreak] [,/deripple] [,/wave2point_correct] [,/float] ; [,/new_avg] [,/subimgx=subimgx] [,subimgy=subimgy] ; [,sllex=sllex] [,slley=slley] [,/nodata] [,/no_calib] ; [,/original] [,/outminsiz] [,/outflatfits] [,/qstop] ; [,/quiet] [,/verbose] [,n_pixel=n_pixel] ; [,run_time=run_time] ; ; trace_prep, filename, image_no of hourly file, index_out, image_out ; trace_prep, infile, dset_arr, index_out, image_out ; trace_prep, index, data, index_out, image_out ; trace_prep, file_list(indir,infile), ss(0:8), index, image ; (with ss=lindgen(total)) ; ; Example for the beginner: ; trace_prep, file_list(indir,infile), [0,2,4,6,8], index, image ,/unspike ,/outtrfits ,outdir=outdir ; ; INPUTS: ; ** There are 2 methods for calling TRACE_PREP ** ; Case 1. input1 - The input TRACE FITS filelist name(s) ; input2 - The dataset number(s) to extract and process ; ; Case 2. input1 - The index structure for each of the input images ; (e.g., index from output of read_trace with ; option /image) ; input2 - The input data array (cube) ; ; ; OUTPUTS (OPTIONAL): ; index_out - The updated index structure of the input images ; image_out - Processed output TRACE images (data cube). ; Default data type is I*2 ; ; OPTIONAL INPUT KEYWORD PARAMETERS: ; normalize - Set to normalize output image to DN per sec ; (TBD) response_norm - Set if want to include response normalization ; to physical units (TBD) ; darkdir - Directory for ~monthly dark current processed fits images, ; default = "concat_dir('$tdb','tdc_darks')", which ; should work at most sites (including vestige and EOF) ; user_dark - If set, read in a user supplied dark current fits image from ; darkdir as the first image with filename given by ; 'user_dark'. This image is assumed to be sized, binned, ; and pixel aligned already. ; no_darksub - Set to not perform dark subtraction; [default: perform it] ; no_flatfield - Set to not perform flat fielding correction (as done ; prior to V3.0); [default: perform it] ; flatdir - Directory for ~quarterly flat field processed fits images, ; default = "concat_dir('$tdb','tdc_darks')", which ; should work at most sites (including vestige and EOF) ; unspike - Set if 1 pass of unspike wanted to remove radiation-belt/ ; cosmic-ray spikes, or set to 1-3 passes of unspike. This ; method uses convolution and thresholding to remove ; spikes ; and may remove small real features. ; [This calls CCK trace_unspike.pro & trace_cleanjpg.pro] ; (TBD) despike - Set if 1 pass of despike, an alternate unspike, wanted to ; (NOT TURNED ON YET!) remove radiation-belt/cosmic-ray spikes, or set to 1-3 ; passes of despike. This method uses a median filter or ; a statistical correction to remove spikes and may remove ; small real features. (TBD) ; [This calls CJS tracedespike.pro & CCK ; trace_cleanjpg.pro] ; destreak - Set if 1 pass of destreak wanted to remove radiation belt/ ; cosmic-ray streaks (particle hits skimming the ; detector). This method uses convolution and thresholding ; for removal and may remove small real features. ; [This calls CCK trace_destreak.pro] ; deripple - Set if 1 pass of readout noise (ripple) remover is wanted. ; This method identifies and zeroes spikes in FFT. ; [This calls CCK trace_knoise.pro] ; wave2point_correct - set to change the pointing tags (xcen,ycen) ONLY, ; based on the wavelength dependent alignment corrections ; contained in trace_wave2point.pro for subsequent use with ; Dominic Zarro's SSW mapping software, e.g., plot_map. ; ***NOTE that these shifts, which are focus dependent, ; are only good at the wavelength focusing positions. In ; recent years the latter have sometimes not been used to ; conserve the focus mechanism. Therefore the researcher ; should separately check image alignments until this ; procedure has been updated. (RWN Nov. 2005)*** ; float - Set if you want to return floating point. [Default is I*2] ; new_avg - Set to recalculate img_avg, img_dev, img_min, img_max for ; image, this is also done upon a call to write_trace to ; save TRACE data cube ; subimgx - Set to x size of subimage for extraction and image_out size; ; to extract a subimage set this and subimgy, sllex, slley ; subimgy - Set to y size of subimage for extraction and image_out size ; sllex - Set to lower left x position of subimage to be extracted ; slley - Set to lower left y position of subimage to be extracted ; nodata - Set if don't want output array (auto set if lt 3 params) ; no_calib - Set to skip calibration but not various output options ; original - Set if input files need more keywords for dark subtraction ; outminsiz - Set if output datacube size to be minimum input image size, ; larger sized images are rebinned to this size after ; being processed; [default is maximum input image size ; with no rebinning] ; outdir - Destination directory for prepped FITS files ; outtrfits - Set if want to output a TRACE binary extension FITS file ; with 1 or more images ; outflatfits - Set if want to output a 2D flat FITS file (1 per image) ; qstop - Set to stop in this procedure for debuging purposes ; quiet - Set for fewer messages, default is loud ; verbose - Set for lots of messages and intermediate data listings; ; suppresses quiet ; ; OPTIONAL OUTPUT KEYWORD PARAMETERS: ; n_pixel - Number of pixels repaired ; run_time - The run time in seconds for TRACE_PREP ; ; COMMON BLOCKS: none. ; ; RESTRICTIONS: ; If the input consists of more than one image, there is no ; check to make sure that all the images are the same size. ; *** The user must be careful to check this. *** ; In fact each image is processed separately with background ; subtractions adjusted for the image size. The output datacube ; is by default sized for the largest image and zero-filled, with ; the images inserted into the cube afterword. There is an ; option to size to the minimum image size. The output datacube ; will be resized to subimage if the subimage extraction option is used. ; ; Note that the usage of the routines called by the keywords unspike, ; despike, destreak, and/or deripple may remove small real features ; in the image. ; ; If the input is an index structure plus a data cube, the index ; is assumed to contain more than the default set of the keywords ; (use '/image' in READ_TRACE; '/all_tags' may result in more then ; 127 keywords, which may cause problems in IDL versions below 5.0) ; ; PROCEDURE: ; Each raw image of the FITS filelist, or of the index and data array, ; are read in one at a time for processing. Missing pixels are ; replaced with the total image average value. Near saturated pixels ; are replaced with values > 4095, such that after dark subtraction they ; are still > 4100, even for the 4x4 summed values. Then the appropriate ; dark pedestal (ADC offset) and the nominally small dark current are ; subtracted from each entire image. Negative data values are left as ; they are. Next the images are flat field corrected. Options are provided ; to call an unspike routine to remove particle and cosmic ray hits of the ; ccd for EUV images only; call an alternate despike program; call a ; destreak procedure to clean, in EUV images only, streaks resulting from ; a cosmic ray traveling through several pixels as it skims through the ; ccd; and call a herringbone readout noise (~few dn ripple) remover for ; cleaning EUV only image backgrounds that have low intensities. Another ; option is to normalize the image for the actual shutter exposure value. ; A subimage can be extracted from each image as an option also. An ; option for updating some of the index parameters for each modified ; image, e.g., 'img_avg', can be selected. ; ; In the future, corrections to the image may also be made for response ; normalization to physical units. Note that no pinhole leakage of foils ; has been observed, so no corrections are needed at the present time. ; The images are returned as an index structure and a datacube, and ; optionally as a TRACE binary extension FITS file with an updated header, ; 1 or more images per file, and a filename of "tfiyyyymmdd_hhmmss", or as ; a 2D flat FITS file (1 image per file) with a filename of ; "tsiyyyymmdd_hhmmss_a#.fits", where # = sai integer (0, 1, or 2). ; ; ; MODIFICATION HISTORY: ;V1.0 Completed on 15-Apr-98 by R. W. Nightingale based on ; eit_prep.pro and sxt_prep.pro ;V1.1 16-Apr-98 (RWN) - Modified so darkdir can be input, if needed ;V1.2 19-Jun-98 (RWN) - Implemented verbose and retrn for TR_DARK_SUB ; - Modified so user_dark file can be supplied to ; TR_DARK_SUB and dark_raw_image & _index get saved ; - Updated procedure description and in-line comments ; - Added image normalization capability ; - Added counting of saturated pixels ;V1.3 6-Jul-98 (RWN) - Modified keyword info. for default 'darkdir' ;V1.4 22-Sep-98 (RWN) - In call to read_trace use /image keyword inplace of ; /all_tags for shorter index list ;V1.5 26-Sep-98 (RWN) - Added option for extraction of a subimage ; - Made the updating of 'img_avg', etc an option ; - Added pointer to dark images at EOF ;V1.6 12-Oct-98 (RWN) - Modify default darkdir keyword documentation ;V1.7 14-Jan-99 (RWN) - Added keyword original to allow more keywords ; - Modify "wave_num" to keyword "wave_len" ;V1.8 6-Mar-99 (SLF) - Add /WAVE2POINT_CORRECT keyword and logic ;V1.9 13-Sep-99 (RWN) - Added more calls to and used UPDATE_HISTORY and ; GET_HISTORY instead of add_tag ; - Added unspike, destreak, despike, deripple, and quiet ; - Updated header and procedure keyword calling sequence ; - Updated /float option and saturated pixel handling ; - Updated to utilize a file/ss map in reading images ;V2.0 11-Feb-00 (RWN) - Updated call to tr_ext_subimg for summed/binned images ;V2.1 14-Mar-00 (RWN) - Correct call to tr_ext_subimg for index2.tbin_ccd ; - Add /quiet pass through to read_trace ;V3.0 23-Aug-02 (RWN) - Add TR_FLAT_SUB interface and UV flat field correction ; - Update library call for mean and standard deviation ;V3.1 17-Sep-02 (RWN) - Include double as possible image type for data_chk ;V3.2 20-Feb-03 (RWN) - Add EUV flat field correction ;V3.2a 4-Nov-05 (RWN) - Header updates, including wave2point note ;- ;----------------------------------------------------------------------------- ; progver = 'V3.2' prognam = 'TRACE_PREP.PRO' t0 = systime(1) t1 = t0 ; Keep track of running time print, 'running ', prognam, ' ', progver, ' with flat field correction' use_fits = 1 wavenum = ['171', '195', '284', '1216', '1550', '1600', '1700', 'WL'] ; loud = 1 - keyword_set(quiet) verbose = keyword_set(verbose) if (verbose eq 1) then loud = 1 no_darksub = keyword_set(no_darksub) no_flatfield = keyword_set(no_flatfield) original = keyword_set(original) if (n_elements(outdir) le 0) then outdir = curdir() yes_data = (not keyword_set(nodata)) and (n_params() gt 3) ; if (datatype(input1) eq 'STC') then begin ; if a structure, need data array of at least 8x8 if (n_elements(input2) gt 63) then use_fits = 0 else begin message, /info, 'No Data supplied, reading FITS file...' ; or read FITS file endelse endif nfiles = n_elements(input1) if (nfiles lt 1) then begin message, /continue, '*** No filename or structure supplied, returning...***' return endif if (use_fits) then begin ; if a FITS file, need image number(s), default = 0 nimages = n_elements(input2) if (nimages eq 0) then begin nimages = 1 input2 = 0 endif ss = input2 read_trace, input1, input2, indices, data, /image, /nodata ; read only the FITS filelist structures ; for sizing the output datacube ; /image must be used to get keywords for dark subtraction ; /all_tags provides too many tags that may run over the IDL limit of 127 for versions below 5.0 if (input2(0) eq -1) then begin ; if reading all images nimages = n_elements(indices) ss = indgen(nimages) endif if (nimages eq 1) then begin ; if 1 image nx0 = gt_tagval(indices,/naxis1) ny0 = gt_tagval(indices,/naxis2) endif else begin ; if more than 1 image if (keyword_set(outminsiz)) then begin ; option for minimum input image size for datacube nx0 = min(gt_tagval(indices,/naxis1)) ny0 = min(gt_tagval(indices,/naxis2)) endif else begin ; default to maximum input image size for datacube nx0 = max(gt_tagval(indices,/naxis1)) ny0 = max(gt_tagval(indices,/naxis2)) endelse endelse print, 'Reading FITS filelist: ',input1,' with image_no(s) = ', input2 if (input2(0) eq -1) then print, ' The total number of images to be read = ',nimages nx2 = nx0 ny2 = ny0 if (keyword_set(subimgx)) then nx2 = subimgx if (keyword_set(subimgy)) then ny2 = subimgy print, 'Output datacube size will be = [', nx2, ',', ny2, ']' endif else begin ; not a FITS file nimages = n_elements(input1) if (nimages eq 1) then begin ; if 1 image nx0 = gt_tagval(input1,/naxis1) ny0 = gt_tagval(input1,/naxis2) endif else begin ; if more than 1 image if (keyword_set(outminsiz)) then begin ; option for minimum input image size for datacube nx0 = min(gt_tagval(input1,/naxis1)) ny0 = min(gt_tagval(input1,/naxis2)) endif else begin ; default to maximum input image size for datacube nx0 = max(gt_tagval(input1,/naxis1)) ny0 = max(gt_tagval(input1,/naxis2)) endelse endelse ss = indgen(nimages) print, 'Read image structure and data for ',nimages,' images' print, 'Output datacube size will be = [', nx0, ',', ny0, ']' endelse indices = 0 ;erase indices as will be input again with data below if (not yes_data) then print, '**No images output to datacube' ; due to /nodata or n_params() lt 4 if (keyword_set(no_calib)) then print,'**No calibation done to images, i.e., no darks subtraction, etc.' if (use_fits) then begin ; make file/ss map to allow input from several hourly files mxf_read_header, input1, pheader, /primary dsmap=mxf_dset_map(pheader,ss) endif for i=0, nimages-1 do begin ; MAIN LOOP: begins ; process 1 image at a time if (use_fits) then begin image_no = dsmap(i).dset file1 = input1(dsmap(i).ifile) if (not loud) then read_trace, file1, image_no, indexi, data, /image, $ original=original, /quiet else $ read_trace, file1, image_no, indexi, data, /image, original=original ; read FITS image and ; at a sublevel rotates image so N is up! index0 = temporary(indexi) nx = gt_tagval(index0,/naxis1) ny = gt_tagval(index0,/naxis2) image0 = temporary(data) ni = n_elements(image0) endif else begin ; else read index structure and data array index0 = input1(i) image_no = i nx = gt_tagval(index0,/naxis1) ny = gt_tagval(index0,/naxis2) image0 = input2(*,*,i) ni = n_elements(image0) if (ni gt nx*ny) then begin image0 = input2(0:nx-1,0:ny-1,i) ni = n_elements(image0) endif endelse if not required_tags(index0,'img_time,amp,sum_ccdx,bin_ccd,tbin_ccd,t_ccd_a,sri_llex,sri_lley,sht_nom', $ missing=missing) then begin box_message,['The following missing tags are required',missing] box_message,'Use /original switch in trace_prep or /image,/original in read_trace' return endif if (ni lt 1) then begin ; if no elements in image print, 'No elements in IMAGE #', i help, /str, index0 ; list the image index message, /continue, 'No elements in IMAGE, skip to next one' endif else begin ; image has elements if (use_fits) then begin if (i eq 0) then begin file2 = file1 print, 'File = ', file2 endif else if (file1 ne file2) then begin file2 = file1 print, 'File = ', file2 endif endif nx2 = nx ny2 = ny if (keyword_set(subimgx)) then nx2 = subimgx if (keyword_set(subimgy)) then ny2 = subimgy wave_numi = gt_tagval(index0,/wave_len) print, 'For image/ ndex',image_no, '/', i, ': Output image size will be = [', nx2, ',', ny2, '], wavelength = ',wave_numi if (keyword_set(verbose)) then $ help, /str, index0 ; list the image index ; CALIBRATION: [default] replace missing and saturated pixels, dark subtract, ; [optional] unspike, destreak, deripple normalization, etc: if (not keyword_set(no_calib)) then begin ; calibrate ; REPLACE MISSING PIXELS WITH IMAGE AVERAGE AND CHECK FOR SATURATED PIXELS ; check .HISTORY (only apply the correction one time) previous = get_history(index0,caller='TRACE_DARK_SUB',found=found) if (not found) then begin mmp = where(image0 eq 0, mcnt) ; find missing pixels if (mcnt ne 0) then begin ; replace if there are any missing nnp = where(image0 gt 0, ncnt) ; find pixels for average if (ncnt ne 0) then begin mean = total(image0(nnp))/ni ; find mean of image if (mean gt 32767. or mean lt -32768.) then begin ; check if I*2 message, /continue, 'Unexpected image mean value outside of I*2: mean = 0.' mean = 0.0 endif endif else message, /info, 'All pixels = 0' ; all pixels = 0 image0(mmp) = fix(mean) ; replace missing pixels with I*2 mean value n_pixel = mcnt ; Update .HISTORY tag tagval0 = string('Replaced ', n_pixel, ' missing pixels with IMAGE_AVERAGE = ', $ mean, format='(a,i5,a,i5)') if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval0 update_history,index0,tagval0 endif ssp = where(image0 ge 4000, scnt) ; find pixels near saturation if (scnt gt 0) then begin image0(ssp) = 4400 ; replace saturated pixels with large value, so that ; 1x1, 2x2 & 4x4 are still >4095 after dark subtraction ; Update .HISTORY tag tagval1 = string('Replaced ', scnt, ' saturated pixels with value > 4100', $ format='(a,i5,a)') if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval1 update_history,index0,tagval1 endif endif ;**** SUBTRACT DARK PEDESTAL (ADC offset) and DARK CURRENT (small) (unless no_darksub set) ; check .HISTORY (only apply the correction one time) previous = get_history(index0,caller='TR_DARK_SUB',found=found) if (not found and not no_darksub) then begin image1 = tr_dark_sub(index0, image0, nx=nx, ny=ny, darkdir=darkdir, $ dark_raw_image=dark_raw_image, dark_raw_index=dark_raw_index, $ version=version, retrn=retrn, user_dark=user_dark, loud=loud, verbose=verbose) ; Update .HISTORY tag tagval2 = 'Subtracted dark pedestal and current ' if (retrn eq 0) then begin if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval2 update_history,index0,tagval2,caller='TR_DARK_SUB',version=version endif else print, '*** NO dark pedestal subtracted' endif else begin if (loud and found) then $ box_message,'TR_DARK_SUB correction already applied, skipping...' image1 = image0 endelse index1 = index0 ;**** Divide image by FLAT FIELD CORRECTION (unless no_flatfield set) ; check .HISTORY (only apply the correction one time) previous = get_history(index1,caller='TR_FLAT_SUB',found=found) if (not found and not no_flatfield) then begin image2 = tr_flat_sub(index1, image1, nx=nx, ny=ny, flatdir=flatdir, $ verbose=verbose, loud=loud, version=version, retrn=retrn, $ ff_1700_image=ff_1700_image, ff_1700_index=ff_1700_index, $ ff_1600_image=ff_1600_image, ff_1600_index=ff_1600_index, $ ff_1550_image=ff_1550_image, ff_1550_index=ff_1550_index, $ ff_1216_image=ff_1216_image, ff_1216_index=ff_1216_index, $ ff_euv_image=ff_euv_image, ff_euv_index=ff_euv_index, $ ff_wl_image=ff_wl_image, ff_wl_index=ff_wl_index) ; Update .HISTORY tag tagval12 = 'Image divided by corrected flat field' if (retrn eq 0) then begin if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval12 update_history,index1,tagval12,caller='TR_FLAT_SUB',version=version endif else print, '*** NO flat field correction applied' endif else begin if (loud and found) then $ box_message,'TR_FLAT_SUB correction already applied, skipping...' image2 = image1 endelse index2 = index1 ; Unspike an EUV image to remove bright bad pixels due to rad-belt/cosmic-ray hits if (index2.wave_len eq '171' or index2.wave_len eq '195' or index2.wave_len eq '284') then begin if (keyword_set(unspike)) then begin if (n_elements(unspike) gt 0) then begin for j = 1, unspike do begin image2 = trace_unspike(temporary(image2), /cleanjpg) endfor tagval3 = string('Unspiked EUV image ',unspike,' time(s) ', format='(a,i3,a)') endif else begin image2 = trace_unspike(temporary(image2), /cleanjpg) tagval3 = 'Unspiked EUV image one time ' endelse ; Update .HISTORY tag if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval3 update_history,index2,tagval3,caller='TRACE_UNSPIKE' endif ; Despike Alternative to EUV unspike if (keyword_set(despike)) then begin ;alternate despike method if (n_elements(unspike) gt 0) then begin for j = 1, unspike do begin image2 = trace_cleanjpg(tracedespike(temporary(image2))) endfor tagval4 = string('Despiked EUV image ',despike,' time(s) ', format='(a,i3,a)') endif else begin image2 = trace_cleanjpg(tracedespike(temporary(image2))) tagval4 = 'Despiked image one time ' endelse ; Update .HISTORY tag if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval4 update_history,index2,tagval4,caller='TRACEDESPIKE' endif ; Destreak an EUV image to remove bright bad streaks due to rad-belt/cosmic-ray hits if (keyword_set(destreak)) then begin image2 = trace_destreak(temporary(image2)) ; Update .HISTORY tag tagval5 = 'Destreaked image ' if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval5 update_history,index2,tagval5,caller='TRACE_DESTREAK' endif ; Remove readout noise (~few dn ripple) from EUV image if (keyword_set(deripple)) then begin image2 = trace_knoise(temporary(image2)) ; Update .HISTORY tag tagval6 = 'Removed readout noise from image ' if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval6 update_history,index2,tagval6,caller='TRACE_KNOISE' endif endif ; Extract a subimage if /subimgx is set to new x axis: ; note: .HISTORY is checked and updated within fuction if (keyword_set(subimgx)) then begin dark = 0 image3 = tr_ext_subimg(index2, image2, index3, nx1=subimgx, dark=dark, $ ny1=subimgy, sri_llex0=sllex, sri_lley0=slley, $ tsum0=index2.tbin_ccd, verbose=verbose) nx0 = subimgx ny0 = subimgy nx = subimgx ny = subimgy tagval7 = string('Subimage extracted with lower left corner at ', $ sllex,slley,' for size ', nx,' x ', ny, ' and tbin =', $ index2.tbin_ccd, format='(a,i5,",",i5,a,i5,a,i5,a,i2)') if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval7 endif else begin index3 = index2 image3 = image2 endelse ; Renormalize to per sec if /normalize is set: ; check .HISTORY (only apply the correction one time) previous = get_history(index3,caller='norm',found=found) if (not found) then begin if (keyword_set(normalize)) then begin exptime = gt_tagval(index3,/sht_mdur) if (exptime gt 0.) then begin image3 = image3 / exptime ; Update .HISTORY tag tagval8 = 'Exposure normalized (per sec)' if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no,': ',tagval8 update_history,index3,tagval8 endif else print, '*** NO normalization because exposure (=< 0.) = ',exptime endif endif else if loud then $ box_message,'Norm correction already applied, skipping...' ; Update image index parameters if /new_avg is set: if (keyword_set(new_avg)) then begin imin = min(image3, max=imax) if (n_elements(image3) eq 1) then begin idev = 0. iavg = image3(0) endif else begin result = moment(image3, /double, sdev=sdev) iavg = float(result(0)) idev = float(sdev) endelse ; Update .HISTORY tag tagval9 = 'Updated img_avg, _dev, _min, _max' index3.img_min = float(imin) index3.img_max = float(imax) index3.img_avg = iavg index3.img_dev = idev if (loud or i eq 0) then $ print, 'HISTORY record updated for image',image_no, $ ': Updated img_avg, _dev, _min, _max' update_history,index3,tagval9 endif endif ; end of calibration ; Return output image as Integer*2 unless /float is set. case data_chk(image3,/type) of 2: if (keyword_set(float)) then image3 = float(temporary(image3)) 4: if (1-keyword_set(float)) then image3 = fix(round(temporary(image3))) 5: if (1-keyword_set(float)) then image3 = fix(round(temporary(image3))) else image3 = float(temporary(image3)) endcase ; Concatenate the headers and images if (i eq 0) then index_out = index3 else begin out_str = index_out(0) index_out = [str_copy_tags(out_str,index_out), str_copy_tags(out_str,index3)] endelse if (i eq 0 and yes_data) then begin if (keyword_set(float)) then type = 4 else type = 2 ; default (=2) is I*2; (=4) is float image_out = make_array(nx0, ny0, nimages, type=type) ; define output datacube endif if (max(image3) ne 0) then begin ; skip blank images if (yes_data) then begin if (nx gt nx0 or ny gt ny0) then begin ; if image size > datacube size image3 = rebin(temporary(image3),nx0,ny0) ; rebin if image integral multiple of datacube print,'For image',image_no,': Rebin image to smaller datacube 2D size' endif image_out(0,0,i) = image3(*,*) ; insert processed image into datacube if (loud or i eq 0) then print,'For image',image_no,': Processed and Inserted into output datacube' endif endif else print, '*** NO output of this zero image for index = ',i endelse ; image has elements systm = systime(1) loop_time = systm - t1 t1 = systm if (loud or i eq 0) then message, /info, string(' 1 image took ', loop_time, ' seconds', $ format='(a,f5.1,a)') endfor ; MAIN LOOP ends if keyword_set(wave2point_correct) then $ index_out=trace_wave2point(index_out) ; apply wave dependent correct ; ; to XCEN/YCEN for latter use with ; ; Dominic Zarro's SSW mapping software. ; Write the optional FITS file if (keyword_set(outtrfits)) then begin write_trace,index_out,image_out,outdir=outdir,prefix='tfi',/loud endif if (keyword_set(outflatfits)) then begin write_trace,index_out,image_out,outdir=outdir,prefix='tsi',/flat_fits,/loud endif run_time = systime(1)-t0 message, /info, string('Total processing time of ', run_time,' seconds for ', $ nimages, ' images' ,format='(a,f8.1,a,i5,a)') if (keyword_set(qstop)) then STOP end