SUBROUTINE DRIZZE C++ C C DRIZZLE.F Version 1.2 Variable Pixel Linear Reconstruction C C Also widely known as "drizzling" C C Earlier versions were known as DRIP. C C This task implements a simple method for transforming one image C onto the pixel grid of another. The method calculates where the C input pixels map on the output grid and performs an optimal weighted C combination of the input and output. Another image records the weighting C at a given pixel. C C It is based largely on ideas of Andy Fruchter and was coded by Richard Hook. C (Contacts: rhook@eso.org & fruchter@stsci.edu) C C This task uses coefficients held in a text file to define C geometrical distortions. C C This code uses the F77/VOS for IRAF i/o. C It may be compiled and linked using a command like: C C xc -p stsdas (-z) drizzle.f drutil.f x_drizzle.x libiraf77.a C C Note that the optional -z switch (which tells IRAF to create a static C executable which doesn't use the IRAF sharable library) may be C necessary to avoid a virtual memory use limitation of around 167MB in C IRAF versions up to 2.10.4. C C History: C C First try, Richard Hook, ST-ECF, October 1995 C Generalised coefficients added, RNH, Nov 95 C Working on V0.2, RNH, Nov 1995 (including "small|large" pixels) C Modified again for continuously variable input pixel area, 29/11/95 C Modified again for X,Y position processing, 1/12/95 C Started work on V0.3 including "update" option, more header C information and WCS, 12/12/95 C V0.3 released, 15th December 1995 C V0.31 modifications following comments from STScI, 27/12/95 C V0.32 modifications, 29/12/95 C V0.33 modifications, STScI, 5th Jan 1996 C V0.34 modifications, try to minimise memory use, 7th Jan 1996 C V0.35 tidyup and a few enhancements, 9th January 1996 C V0.36 correct rotation angle in output and also case C where outdat=olddat, 8th February 1996 C V0.37 correct silly rotation angle multiplication error C 29th February 1996 C V0.40 first pre-release version, tidied up and corrected WCS C bug with CTYPE header items. 25th March 1996 C V0.5 started to incorporate the "boxer" code for linear C combination by precise area weighting rather than by the C drizzling approximation. The "boxer" code was originally C written by Bill Sparks, STScI, 19th April 1996 C V0.9 many simpifications and improvements for a pre-release C version. This removes drizzling and the XY option as well C as many options which weren't used. 16th May 1996 C V0.91 corrected typo when shunits=output. 21st May 1996 C V0.92 changed parameter names and removed "olddata" option. 22nd May 1996 C V0.93 changed to separate coeffients files for each coefficients C set. 22nd May 1996 C V0.94 changed access to coefficients file to use IRAF filenames via C (hidden) F77 VOS calls and also put in a scale for the weights C image and changed the names. ST-ECF, 30th May 1996 C V0.95 bug fix in error message within SGAREA, format statement syntax C error. ST-ECF, 20th June 1996 C V0.96 name change to DRIZZLE. ST-ECF, 5th July 1996 C V0.97 added exposure time keyword handling. ST-ECF, 9th July 1996 C V0.98 added ability NOT to store output counts image, 16th July 1996 C V0.99 added flag for scaling by exposure time, 16th July 1996 C V0.991 added "fillval" parameter, option for wt_scl=exptime, 29th July 1996 C V0.992 slight change to fillval interpretation, 31st July 1996 C V0.993 small bugfix when reading coefficients, 28th August 1996 C C V1.0 release. 17th September 1996. C C V1.01 - some bug fixes after first release, 19th September 1996 C V1.02 - same as V1.01 except permissions of tar file altered, 3rd October 1996 C V1.03 - small change to WCS logic to avoid crashes, 25th October 1996 C V1.08 - collection of several small changes, 5th December 1996 C (including a change in the definition of the centres and the C addition of the out_typ parameter)) C V1.09 - correct scaling of offsets bug, 16th December 1996 C V1.091 - introduce Jacobian for weights, 10th January 1997 C V1.092 - correct header update bug, 21st January 1997 C V1.093 - put utility routines in separate file (drutil.f), 4th February 1997 C V1.094 - change disk access to minimise memory use, 10th February 1997 C V1.095 - introduced automatic coeffs/lambda from header, 7th March 1997 C V1.096 - added tidier memory management in error cases, 10th March 1997 C also better file closure on error trapping. C V1.097 - bugfix to V1.096, 12th March 1997 C V1.098 - pre-release version incorporating EXPSQ weighting, modified C error messages and a maximum of 999 images, 22nd April 1997 C C V1.1 - second public release. 12th May 1997, STScI C C V1.2 - first release within STSDAS dither package, this version uses C the in_un and out_un system for defining whether input/outputs C are in count or counts/s. The old parameter exp_sc has been C removed. 20th February 1998. C C-- IMPLICIT NONE C Iraf global storage (real in this program) REAL MEMR(1) !I COMMON /MEM/MEMR !I C Local variables INTEGER ISTAT,I,LUN INTEGER NDIMS,DNX,DNY,ONX,ONY,DDIMS(7),ODIMS(7),DATTYP INTEGER IDD,IDW,IDND,IDNC,IC INTEGER PDATA,PWEI,PNDAT,PNCOU REAL XCO(10),YCO(10),XSH,YSH,LAM,WTSCL REAL SCALE,ROT,PFRACT REAL EXPIN,EXPOUT,EXPCUR,FILVAL REAL LASTPF,LASTSC,LASTOX,LASTOY,OX,OY DOUBLE PRECISION WCS(8) CHARACTER*80 DATA,WEIGHT,OUTDAT,OUTCOU,CHARS CHARACTER*80 COEFFS,SHFTUN,WTSTR,FILSTR CHARACTER*45 VERS CHARACTER*8 CTYPE1,CTYPE2,LASTUN CHARACTER*8 EXPKEY,SHFTFR,INUN,OUTUN,ALIGN C Logical flags LOGICAL VERBOSE LOGICAL NOCO LOGICAL USEWEI LOGICAL UPDATE LOGICAL GOTWCS LOGICAL WRIWEI LOGICAL FILL LOGICAL ROTFIR LOGICAL INCPS LOGICAL OUTCPS LOGICAL CURCPS LOGICAL HCOEFF LOGICAL ADATA LOGICAL AWEI LOGICAL ANDAT LOGICAL ANCOU LOGICAL ODD LOGICAL ODW LOGICAL ODND LOGICAL ODNC C-- Start of executable code VERS='DRIZZLE Version 1.2 (February 1998)' C First announce the version CALL UMSPUT('+ '//VERS,1,0,ISTAT) VERBOSE=.TRUE. C Also get fraction of area of input pixels to be used CALL UCLGSR('pixfrac',PFRACT,ISTAT) C Check for small values IF(PFRACT.LT.0.001) PFRACT=0.001 C Get the filling value for undefined pixels CALL UCLGST('fillval',FILSTR,ISTAT) C Get the switch value for whether to read in CPS or C counts (latter is the default) CALL UCLGST('in_un',INUN,ISTAT) IF(INUN.EQ.'cps') THEN INCPS=.TRUE. ELSE INCPS=.FALSE. ENDIF C Get the switch value for whether to write out in CPS or C counts (latter is the default) CALL UCLGST('out_un',OUTUN,ISTAT) IF(OUTUN.EQ.'cps') THEN OUTCPS=.TRUE. ELSE OUTCPS=.FALSE. ENDIF C Get parameter specifying which frames shifts are applied CALL UCLGST('shft_fr',SHFTFR,ISTAT) IF(SHFTFR.EQ.'input') THEN ROTFIR=.FALSE. ELSE ROTFIR=.TRUE. ENDIF C Check for meaningful values IF(FILSTR.NE.'INDEF'.AND.FILSTR.NE.'indef') THEN READ(FILSTR,*,IOSTAT=ISTAT) FILVAL IF(ISTAT.NE.0) THEN CALL UMSPUT('! Invalid filling value specified', : 1,0,ISTAT) GO TO 99 ENDIF FILL=.TRUE. ELSE FILL=.FALSE. ENDIF C Get the origin of the image coordinates - in other words C whether we treat the center or the corner of the central C pixel as the reference point CALL UCLGST('align',ALIGN,ISTAT) C Get the name of the file containing the coefficients CALL UCLGST('coeffs',COEFFS,ISTAT) C Get the scale and rotation for the final image CALL UCLGSR('scale',SCALE,ISTAT) CALL UCLGSR('rot',ROT,ISTAT) C Check for invalid scale IF(SCALE.EQ.0.0) THEN CALL UMSPUT('! Invalid scale',1,0,ISTAT) GO TO 99 ENDIF C Convert the rot to radians ROT=ROT*3.1415926/180.0 C Get the X,Y final shift CALL UCLGSR('xsh',XSH,ISTAT) CALL UCLGSR('ysh',YSH,ISTAT) C Also get the binary flag for whether units are input or C output pixels CALL UCLGST('shft_un',SHFTUN,ISTAT) IF(SHFTUN.EQ.'output') THEN XSH=XSH*SCALE YSH=YSH*SCALE ENDIF C Note that the shifts we now have are INPUT ones C Get the keyword used for exposure time CALL UCLGST('expkey',EXPKEY,ISTAT) C Check for compatible requests IF(EXPKEY.EQ.' ') THEN CALL UMSPUT('! No exposure time keyword has been given', : 1,0,ISTAT) GO TO 99 ENDIF C Get the name of the input data image CALL UCLGST('data',DATA,ISTAT) C Try to access the input data image CALL UIMOPN(DATA,1,IDD,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to open data image', : 1,0,ISTAT) ODD=.FALSE. GO TO 99 ELSE IF(VERBOSE) CALL UMSPUT( : '-Opening input data image: '//DATA, : 1,0,ISTAT) ODD=.TRUE. ENDIF C Get the array size and shape CALL UIMGID(IDD,DATTYP,NDIMS,DDIMS,ISTAT) IF(NDIMS.NE.2) THEN CALL UMSPUT('! Input image is not two dimensional', : 1,0,ISTAT) GO TO 99 ELSE DNX=DDIMS(1) DNY=DDIMS(2) ENDIF C Interpret the coefficients file name and act appropriately NOCO=.FALSE. IF(COEFFS.EQ.' ') NOCO=.TRUE. C Now check for the special string "header" IF(COEFFS.EQ.'header') THEN HCOEFF=.TRUE. C and try to read these from the input header CALL GTCOEF(IDD,COEFFS,LAM,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to get necessary header information', : 1,0,ISTAT) GO TO 99 ENDIF ELSE HCOEFF=.FALSE. ENDIF C Set default values IF(NOCO) THEN DO I=1,10 XCO(I)=0.0 YCO(I)=0.0 ENDDO ELSE C Open the coefficients file CALL UFOPEN(COEFFS,1,LUN,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to open coefficients file',1,0, : ISTAT) GO TO 99 ELSE IF(VERBOSE) CALL UMSPUT( : '-Opening coefficients file: '//COEFFS, : 1,0,ISTAT) ENDIF C Get the wavelength (nm). This is only used for "trauger" C style coefficients sets, if we are using the header we C will already have this IF(.NOT.HCOEFF) CALL UCLGSR('lambda',LAM,ISTAT) C Read in the coefficients from the file CALL GETCO(LUN,LAM,XCO,YCO,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to read coefficients', : 1,0,ISTAT) CALL UFCLOS(LUN,ISTAT) GO TO 99 ENDIF C Close the coefficients file CALL UFCLOS(LUN,ISTAT) ENDIF C Get the input exposure time from the header CALL UHDGSR(IDD,EXPKEY,EXPIN,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Error, unable to read exposure time keyword value', : 1,0,ISTAT) CALL UMSPUT( : '! Please add one to the image header and try again.', : 1,0,ISTAT) GO TO 99 ENDIF C Try to also get a "weighting" image and a scale factor for it CALL UCLGST('in_mask',WEIGHT,ISTAT) CALL UCLGST('wt_scl',WTSTR,ISTAT) C The scale is read as a string as it might be "exptime" C in which case we use the exposure time - if we have it C We also can have weighting by exposure time squared (eg, C for read-noise dominated images) IF(WTSTR.EQ.'exptime' .OR. WTSTR.EQ.'EXPTIME') THEN WTSCL=EXPIN ELSE IF(WTSTR.EQ.'expsq' .OR. WTSTR.EQ.'EXPSQ') THEN WTSCL=EXPIN**2 ELSE READ(WTSTR,*,IOSTAT=ISTAT) WTSCL IF(ISTAT.NE.0) THEN CALL UMSPUT('! Weighting scale is not a number', : 1,0,ISTAT) GO TO 99 ENDIF ENDIF IF(WEIGHT.EQ.' ') THEN USEWEI=.FALSE. ODW=.FALSE. ELSE CALL UIMOPN(WEIGHT,1,IDW,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Unable to open mask image', : 1,0,ISTAT) ODW=.FALSE. GO TO 99 ELSE IF(VERBOSE) CALL UMSPUT( : '-Opening mask image: '//WEIGHT, : 1,0,ISTAT) ODW=.TRUE. C Get the array size and shape of the mask image - must C be the same as the data CALL UIMGID(IDW,DATTYP,NDIMS,DDIMS,ISTAT) IF(NDIMS.NE.2 .OR. DDIMS(1).NE.DNX .OR. : DDIMS(2).NE.DNY) THEN CALL UMSPUT('! Mask image is not the same shape '// : 'and size as data',1,0,ISTAT) GO TO 99 ENDIF USEWEI=.TRUE. ENDIF ENDIF C We get the name of the output images here CALL UCLGST('outdata',OUTDAT,ISTAT) CALL UCLGST('outweig',OUTCOU,ISTAT) C If the output weighting image is null we just don't C create anything and don't write anything out IF(OUTCOU.EQ.' ') THEN WRIWEI=.FALSE. ELSE WRIWEI=.TRUE. ENDIF C It may be that the "new" image currently exists and we C want to update it by drizzling the new data on top IF(OUTDAT.EQ.' ') THEN CALL UMSPUT( : '! No output data array name specified', : 1,0,ISTAT) GO TO 99 ENDIF C We try to open it with update access CALL UIMOPN(OUTDAT,2,IDND,ISTAT) IF(ISTAT.NE.0) THEN ODND=.FALSE. UPDATE=.FALSE. ELSE CALL UMSPUT( : '-Opening extant data image: '//OUTDAT, : 1,0,ISTAT) ODND=.TRUE. C Now we try to open the corresponding weighting image CALL UIMOPN(OUTCOU,2,IDNC,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to open extant weighting image', : 1,0,ISTAT) ODNC=.FALSE. GO TO 99 ENDIF CALL UMSPUT( : '-Opening extant weighting image: '//OUTCOU, : 1,0,ISTAT) ODNC=.TRUE. C We are doing well and probably can use the existing C images however we need to do a little more... C Get the array size and shape of old data frame CALL UIMGID(IDND,DATTYP,NDIMS,ODIMS,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! A problem was encountered accessing existing image: '// : OUTDAT,1,0,ISTAT) CALL UMSPUT( : '! Please delete or recreate it and try again', : 1,0,ISTAT) GO TO 99 ENDIF IF(NDIMS.NE.2) THEN CALL UMSPUT('! Extant image is not two dimensional', : 1,0,ISTAT) GO TO 99 ELSE ONX=ODIMS(1) ONY=ODIMS(2) ENDIF C Also check that the extant weighting image is the same size CALL UIMGID(IDNC,DATTYP,NDIMS,ODIMS,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! A problem was encountered accessing existing image: '// : OUTDAT,1,0,ISTAT) CALL UMSPUT( : '! Please delete or recreate it and try again', : 1,0,ISTAT) GO TO 99 ENDIF IF(NDIMS.NE.2) THEN CALL UMSPUT( : '! Current weighting image is not two dimensional', : 1,0,ISTAT) GO TO 99 ELSE IF(ODIMS(1).NE.ONX .OR. ODIMS(2).NE.ONY) THEN CALL UMSPUT('! Count image not the same size as data', : 1,0,ISTAT) GO TO 99 ENDIF C We can now set the update flag UPDATE=.TRUE. WRIWEI=.TRUE. C Get the current exposure time from the header CALL UHDGSR(IDND,EXPKEY,EXPCUR,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Warning, failed to get exposure time keyword value', : 1,0,ISTAT) GO TO 99 ENDIF C We now get some of the previous drizzle information from C the header of the current image and check for funny C combinations CALL GOLDH(IDND,IC,LASTPF,LASTSC, : LASTUN,LASTOX,LASTOY,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Warning, unable to read old drizzle header keywords', : 1,0,ISTAT) CALL UMSPUT('! assuming current image is in counts', : 1,0,ISTAT) CURCPS=.FALSE. ELSE WRITE(CHARS, : '(''-This is image # '',I3,'' to be drizzled onto this output''' : //')') IC+1 CALL UMSPUT(CHARS,1,0,ISTAT) IF(LASTUN.EQ.'counts' .OR. LASTUN.EQ.'COUNTS') THEN CURCPS=.FALSE. ELSE CURCPS=.TRUE. ENDIF IF(LASTPF.NE.PFRACT) THEN CALL UMSPUT( : '! Warning, current pixfrac does not match previous one', : 1,0,ISTAT) ENDIF IF(LASTSC.NE.SCALE) THEN CALL UMSPUT( : '! Warning, current scale does not match previous one', : 1,0,ISTAT) ENDIF IF(ALIGN.EQ.'corner') THEN OX=REAL(ONX/2)+0.5 OY=REAL(ONY/2)+0.5 ELSE OX=REAL(ONX/2)+1.0 OY=REAL(ONY/2)+1.0 ENDIF IF(LASTOX.NE.OX .OR. LASTOY.NE.OY) THEN CALL UMSPUT( : '! Warning, current reference pixel alignment '// : 'does not match previous one', : 1,0,ISTAT) ENDIF ENDIF ENDIF C If there is no extant output image C we will have to create output data and weighting images anyway C so we need to find out their dimensions IF(.NOT.UPDATE) THEN C Read the dimensions for the output image from the CL, if it C is to be newly created CALL UCLGSI('outnx',ONX,ISTAT) CALL UCLGSI('outny',ONY,ISTAT) ODIMS(1)=ONX ODIMS(2)=ONY C If we are not updating we need to create the output images CALL UIMCRE(OUTDAT,6,2,ODIMS,IDND,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to create output data image', : 1,0,ISTAT) ODND=.FALSE. GO TO 99 ELSE ODND=.TRUE. ENDIF C If we are going to store the counts image we must create it IF(WRIWEI) THEN CALL UIMCRE(OUTCOU,6,2,ODIMS,IDNC,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to create output weighting image', : 1,0,ISTAT) ODNC=.FALSE. GO TO 99 ELSE ODNC=.TRUE. ENDIF ELSE CALL UMSPUT( : '! Warning, output weighting image will not be saved', : 1,0,ISTAT) ENDIF ENDIF C If this is a new output then we need to copy the C old header over IF(.NOT.UPDATE) THEN CALL UHDCPY(IDD,IDND,ISTAT) IF(WRIWEI) CALL UHDCPY(IDD,IDNC,ISTAT) ENDIF C Update the exposure times in the outputs IF(UPDATE) THEN EXPOUT=EXPIN+EXPCUR ELSE EXPOUT=EXPIN ENDIF CALL UHDPSR(IDND,EXPKEY,EXPOUT,ISTAT) IF(WRIWEI) CALL UHDPSR(IDNC,EXPKEY,EXPOUT,ISTAT) C Now we allocate space for all the arrays we will need, note that they C are all real C C First one line of the data array CALL UDMGET(DNX,6,PDATA,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to allocate memory for data array line', : 1,0,ISTAT) ADATA=.FALSE. GO TO 99 ELSE ADATA=.TRUE. ENDIF C Mask array - just one line at a time CALL UDMGET(DNX,6,PWEI,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to allocate memory for weights array', : 1,0,ISTAT) AWEI=.FALSE. GO TO 99 ELSE AWEI=.TRUE. ENDIF C Weighting output/current array CALL UDMGET(ONX*ONY,6,PNCOU,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to allocate memory for output weighting array', : 1,0,ISTAT) ANCOU=.FALSE. GO TO 99 ELSE ANCOU=.TRUE. ENDIF C Output/current data array CALL UDMGET(ONX*ONY,6,PNDAT,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Unable to allocate memory for output data array', : 1,0,ISTAT) ANDAT=.FALSE. GO TO 99 ELSE ANDAT=.TRUE. ENDIF C Also try to get the header items for the WCS, note that C this is not the full WCS but just 8 header items commonly used C for HST images plus CTYPE1 & CTYPE2 GOTWCS=.FALSE. IF(.NOT.UPDATE) THEN CALL GETWCS(IDD,WCS,CTYPE1,CTYPE2,ISTAT) IF(ISTAT.EQ.0) THEN GOTWCS=.TRUE. ELSE CALL UMSPUT( : '! Warning, unable to read WCS information from header', : 1,0,ISTAT) ENDIF ENDIF C Here we set the WCS matrix to something sensible anyway IF(.NOT.GOTWCS) THEN CTYPE1='PIXELS' CTYPE2='PIXELS' DO I=1,8 WCS(I)=0.0D0 ENDDO ENDIF C Optional weights array IF(.NOT.USEWEI) CALL SETIM(MEMR(PWEI),DNX,1,1.0) C If we are updating we read in the extant "new" images IF(UPDATE) THEN DO I=1,ONY CALL UIGL2R(IDND,I,MEMR(PNDAT+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Failed to read current data image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO C Scale by exposure time IF(.NOT.CURCPS) THEN CALL MULC(MEMR(PNDAT),ONX,ONY,1.0/EXPCUR) ENDIF DO I=1,ONY CALL UIGL2R(IDNC,I,MEMR(PNCOU+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT( : '! Failed to read current weighting image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO ELSE C The final case is where the output is empty. C Before we start drizzling we need to drain out the output arrays C if there are no old arrays CALL SETIM(MEMR(PNDAT),ONX,ONY,0.0) CALL SETIM(MEMR(PNCOU),ONX,ONY,0.0) ENDIF C Now do the actual combination using BOXER CALL DOBOX(MEMR(PDATA),MEMR(PWEI), : MEMR(PNDAT),MEMR(PNCOU),DNX,DNY, : ONX,ONY,XCO,YCO,WTSCL,ALIGN,INCPS,EXPIN, : PFRACT,SCALE,ROT,XSH,YSH,WCS,ROTFIR, : UPDATE,IDD,USEWEI,IDW,ISTAT) C Check the exit status IF(ISTAT.NE.0) GO TO 99 C Scale the output if not CPS output IF(.NOT.OUTCPS) CALL MULC(MEMR(PNDAT),ONX,ONY,EXPOUT) C Put in the fill values (if defined) IF(FILL) CALL PUTFIL(MEMR(PNDAT),MEMR(PNCOU),ONX,ONY,FILVAL) C Write out the two new images IF(VERBOSE) : CALL UMSPUT('-Writing output drizzled image: '// : OUTDAT,1,0,ISTAT) DO I=1,ONY CALL UIPL2R(IDND,I,MEMR(PNDAT+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to write output data image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO C Update the header CALL UHEAD(IDND,VERS,DATA,WEIGHT,WTSCL,OUTDAT, : OUTCOU,COEFFS,SHFTUN,PFRACT, : LAM,SCALE,ROT,SHFTFR,ALIGN, : EXPKEY,FILSTR,INUN,OUTUN, : XSH,YSH,DNX,DNY,ONX,ONY) C If we are not updating (ie, this is the first drizzle C onto this output) then we modify the WCS IF((.NOT.UPDATE).AND.GOTWCS) THEN CALL UWCS(IDND,WCS,CTYPE1,CTYPE2) ENDIF C Now write out the weights array, if there is one IF(WRIWEI) THEN IF(VERBOSE) : CALL UMSPUT('-Writing output weighting image: '// : OUTCOU,1,0,ISTAT) DO I=1,ONY CALL UIPL2R(IDNC,I,MEMR(PNCOU+(I-1)*ONX),ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to write output weight image', : 1,0,ISTAT) GO TO 99 ENDIF ENDDO C Update the header CALL UHEAD(IDNC,VERS,DATA,WEIGHT,WTSCL,OUTDAT, : OUTCOU,COEFFS,SHFTUN,PFRACT, : LAM,SCALE,ROT,SHFTFR,ALIGN, : EXPKEY,FILSTR,INUN,OUTUN, : XSH,YSH,DNX,DNY,ONX,ONY) IF((.NOT.UPDATE).AND.GOTWCS) THEN CALL UWCS(IDNC,WCS,CTYPE1,CTYPE2) ENDIF ENDIF 99 CONTINUE C Free up all allocated memory arrays IF(ADATA) CALL UDMFRE(PDATA,6,ISTAT) IF(AWEI) CALL UDMFRE(PWEI,6,ISTAT) IF(ANCOU) CALL UDMFRE(PNCOU,6,ISTAT) IF(ANDAT) CALL UDMFRE(PNDAT,6,ISTAT) C Close all the image files which may have been open IF(ODD) CALL UIMCLO(IDD,ISTAT) IF(ODW) CALL UIMCLO(IDW,ISTAT) IF(ODND) CALL UIMCLO(IDND,ISTAT) IF(ODNC) CALL UIMCLO(IDNC,ISTAT) C End of main DRIZZLE module END SUBROUTINE UHEAD(ID,VERS,DATA,WEIGHT,WTSCL,OUTDAT, : OUTCOU,COEFFS,SHFTUN,PFRACT, : LAM,SCALE,ROT,SHFTFR,ALIGN,EXPKEY, : FILSTR,INUN,OUTUN, : XSH,YSH,DNX,DNY,ONX,ONY) C C Update the header of the output image with all the C parameters of the current DRIZZLE run. C IMPLICIT NONE INTEGER ID,DNX,DNY,ONX,ONY,ISTAT,I CHARACTER*80 DATA,WEIGHT CHARACTER*80 OUTDAT,OUTCOU,COEFFS,SHFTUN CHARACTER*80 BUFFER,FILSTR CHARACTER*40 VERS CHARACTER*8 EXPKEY,INUN,OUTUN,ALIGN CHARACTER*4 STEM CHARACTER*8 SHFTFR REAL PFRACT,LAM,SCALE,ROT,WTSCL REAL DROT,XSH,YSH,XSHI,YSHI,OFF STEM='D000' C First find out if there are old values C This is very ugly DO I=1,999 WRITE(STEM(2:4),'(I3)') I IF(STEM(2:2).EQ.' ') STEM(2:2)='0' IF(STEM(3:3).EQ.' ') STEM(3:3)='0' CALL UHDGST(ID,STEM//'VER',BUFFER,ISTAT) C Check for a status code indicating that this header C item doesn't exist IF(ISTAT.EQ.40) GO TO 111 ENDDO 111 CONTINUE CALL UHDAST(ID,STEM//'VER',VERS, : 'Drizzle, task version',0,ISTAT) C Now we continue to add the other items using the same C "stem" CALL UHDAST(ID,STEM//'DATA',DATA, : 'Drizzle, input data image',0,ISTAT) CALL UHDAST(ID,STEM//'OUDA',OUTDAT, : 'Drizzle, output data image',0,ISTAT) CALL UHDAST(ID,STEM//'OUWE',OUTCOU, : 'Drizzle, output weighting image',0,ISTAT) CALL UHDAST(ID,STEM//'MASK',WEIGHT, : 'Drizzle, input data quality mask image', : 0,ISTAT) CALL UHDASR(ID,STEM//'WTSC',WTSCL, : 'Drizzle, weighting factor for input image',0,ISTAT) CALL UHDASR(ID,STEM//'PIXF',PFRACT, : 'Drizzle, linear size of drop',0,ISTAT) CALL UHDASR(ID,STEM//'SCAL',SCALE, : 'Drizzle, scale (pixel size) of output image',0,ISTAT) CALL UHDAST(ID,STEM//'COEF',COEFFS, : 'Drizzle, coefficients file name ',0,ISTAT) CALL UHDASR(ID,STEM//'LAM',LAM, : 'Drizzle, wavelength applied for transformation (nm)', : 0,ISTAT) C Convert the rotation angle back to degrees DROT=ROT/3.1415926535897*180.0 CALL UHDASR(ID,STEM//'ROT',DROT, : 'Drizzle, rotation angle, degrees anticlockwise',0,ISTAT) C Check the SCALE and units C The units are INPUT pixels on entry to this routine IF(SHFTUN.EQ.'output') THEN XSHI=XSH/SCALE YSHI=YSH/SCALE ELSE XSHI=XSH YSHI=YSH ENDIF CALL UHDASR(ID,STEM//'XSH',XSHI, : 'Drizzle, X shift applied',0,ISTAT) CALL UHDASR(ID,STEM//'YSH',YSHI, : 'Drizzle, Y shift applied',0,ISTAT) CALL UHDAST(ID,STEM//'SFTU',SHFTUN, : 'Drizzle, units used for shifts (output or input pixels)', : 0,ISTAT) CALL UHDAST(ID,STEM//'SFTF',SHFTFR, : 'Drizzle, frame in which shifts were applied', : 0,ISTAT) CALL UHDASI(ID,STEM//'OUNX',ONX, : 'Drizzle, X size of output image (if newly created)',0,ISTAT) CALL UHDASI(ID,STEM//'OUNY',ONY, : 'Drizzle, Y size of output image (if newly created)',0,ISTAT) CALL UHDAST(ID,STEM//'EXKY',EXPKEY, : 'Drizzle, exposure keyword name in input image',0,ISTAT) CALL UHDAST(ID,STEM//'INUN',INUN, : 'Drizzle, units of input image - counts or cps',0,ISTAT) CALL UHDAST(ID,STEM//'OUUN',OUTUN, : 'Drizzle, units of output image - counts or cps',0,ISTAT) CALL UHDAST(ID,STEM//'FVAL',FILSTR, : 'Drizzle, fill value for zero weight output pixels', : 0,ISTAT) IF(ALIGN.EQ.'corner') THEN OFF=0.5 ELSE OFF=1.0 ENDIF CALL UHDASR(ID,STEM//'INXC',FLOAT(DNX/2)+OFF, : 'Drizzle, reference center of input image (X)',0,ISTAT) CALL UHDASR(ID,STEM//'INYC',FLOAT(DNY/2)+OFF, : 'Drizzle, reference center of input image (Y)',0,ISTAT) CALL UHDASR(ID,STEM//'OUXC',FLOAT(ONX/2)+OFF, : 'Drizzle, reference center of output image (X)',0,ISTAT) CALL UHDASR(ID,STEM//'OUYC',FLOAT(ONY/2)+OFF, : 'Drizzle, reference center of output image (Y)',0,ISTAT) RETURN END SUBROUTINE DOBOX(DATA,WEI,NDAT,NCOU,DNX,DNY, : ONX,ONY,XCO,YCO,WTSCL,ALIGN,INCPS,EXPIN, : PFRACT,SCALE,ROT,XSH,YSH,WCS,ROTFIR, : UPDATE,IDD,USEWEI,IDW,ISTAT) C C This module does the actual mapping of input flux to output images using C "boxer", a code written by Bill Sparks for FOC geometric C distortion correction, rather than the "drizzling" approximation. C C This works by calculating the positions of the four corners of C a quadrilateral on the output grid corresponding to the corners C of the input pixel and then working out exactly how much of each C pixel in the output is covered, or not. C IMPLICIT NONE INTEGER DNX,DNY,ONX,ONY,IDD,IDW INTEGER I,J,II,JJ,NMISS,ISTAT,NHIT,IS REAL DATA(DNX),WEI(DNX),W,EXPIN REAL NDAT(ONX,ONY),NCOU(ONX,ONY) REAL XCO(10),YCO(10) REAL A,B,C,D DOUBLE PRECISION WCS(8) REAL X,Y,XOUT(4),YOUT(4),VC,XSH,YSH,DOVER,DH REAL EVALCU,SINTH,COSTH,ROT,XF,YF,XOFF,YOFF,SCALE REAL XS,XC,YS,YC,XT,YT,PFRACT,S2,XNL,YNL,XP,YP REAL XNLD,YNLD,DET REAL XCEN,YCEN,XCORN(4),YCORN(4),AC,WTSCL,JACO CHARACTER*80 CHARS CHARACTER*8 ALIGN LOGICAL NONLIN,ROTFIR,UPDATE,USEWEI,INCPS C-- C Some initial settings - note that the reference pixel C position is determined by the value of ALIGN NMISS=0 IF(ALIGN.EQ.'corner') THEN XCEN=FLOAT(DNX/2)+0.5 YCEN=FLOAT(DNY/2)+0.5 ELSE XCEN=FLOAT(DNX/2)+1.0 YCEN=FLOAT(DNY/2)+1.0 ENDIF C If all the non-linear coefficients are zero we can C ignore that part of things NONLIN=.FALSE. DO I=1,10 IF(XCO(I).NE.0.0 .OR. YCO(I).NE.0.0) NONLIN=.TRUE. ENDDO C Calculate the position of the first whole pixel Y=-YCEN C Calculate some numbers to simplify things later SINTH=SIN(ROT) COSTH=COS(ROT) C The shifts are in units of INPUT pixels. XOFF=XSH/SCALE YOFF=YSH/SCALE XF=1.0/SCALE YF=1.0/SCALE C Some more economies S2=1.0/(XF*YF) XS=XF*SINTH XC=XF*COSTH YS=YF*SINTH YC=YF*COSTH IF(ALIGN.EQ.'corner') THEN XP=REAL(ONX/2)+0.5 YP=REAL(ONY/2)+0.5 ELSE XP=REAL(ONX/2)+1.0 YP=REAL(ONY/2)+1.0 ENDIF XT=XOFF+XP YT=YOFF+YP DH=0.5*PFRACT AC=1.0/(PFRACT*PFRACT) C Outer loop over input image pixels (X,Y) DO J=1,DNY Y=Y+1.0 C Read in one row of data CALL UIGL2R(IDD,J,DATA,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to read data image', : 1,0,ISTAT) RETURN ENDIF C If the input image is not in CPS we need to divide by C the exposure IF(.NOT.INCPS) CALL MULC(DATA,DNX,1,1.0/EXPIN) C If there is one, read in a line of the mask image IF(USEWEI) THEN CALL UIGL2R(IDW,J,WEI,ISTAT) IF(ISTAT.NE.0) THEN CALL UMSPUT('! Failed to read mask image', : 1,0,IS) RETURN ENDIF ENDIF X=-XCEN DO I=1,DNX X=X+1.0 C Calculate where this pixel lands up C The angles are assumed to be in radians at this point IF(NONLIN) THEN XCORN(1)=EVALCU(X-DH,Y+DH,XCO) XCORN(2)=EVALCU(X+DH,Y+DH,XCO) XCORN(3)=EVALCU(X+DH,Y-DH,XCO) XCORN(4)=EVALCU(X-DH,Y-DH,XCO) YCORN(1)=EVALCU(X-DH,Y+DH,YCO) YCORN(2)=EVALCU(X+DH,Y+DH,YCO) YCORN(3)=EVALCU(X+DH,Y-DH,YCO) YCORN(4)=EVALCU(X-DH,Y-DH,YCO) ELSE XCORN(1)=X-DH XCORN(2)=X+DH XCORN(3)=X+DH XCORN(4)=X-DH YCORN(1)=Y+DH YCORN(2)=Y+DH YCORN(3)=Y-DH YCORN(4)=Y-DH ENDIF C Apply the linear transform C There are two ways this can be done - shift then C rotate or rotate then shift, the latter was the C method used in all versions of drizzle before 1.04 C when an option was introduced. IF(ROTFIR) THEN DO II=1,4 XOUT(II)=XC*XCORN(II)-YS*YCORN(II)+XT YOUT(II)=XS*XCORN(II)+YC*YCORN(II)+YT ENDDO ELSE DO II=1,4 XOUT(II)=XC*(XCORN(II)+XSH)-YS*(YCORN(II)+YSH)+XP YOUT(II)=XS*(XCORN(II)+XSH)+YC*(YCORN(II)+YSH)+YP ENDDO ENDIF C Work out the area of the quadrilateral on the output grid C Note that this expression expects the points to be in C clockwise order JACO=0.5*((XOUT(2)-XOUT(4))*(YOUT(1)-YOUT(3)) - : (XOUT(1)-XOUT(3))*(YOUT(2)-YOUT(4))) NHIT=0 C Loop over output pixels which could be affected DO JJ=NINT(MIN(YOUT(1),YOUT(2),YOUT(3),YOUT(4))), : NINT(MAX(YOUT(1),YOUT(2),YOUT(3),YOUT(4))) DO II=NINT(MIN(XOUT(1),XOUT(2),XOUT(3),XOUT(4))), : NINT(MAX(XOUT(1),XOUT(2),XOUT(3),XOUT(4))) C Check it is on the output image IF(II.GE.1 .AND. II.LE.ONX .AND. : JJ.GE.1 .AND. JJ.LE.ONY) THEN C Call boxer to calculate overlap CALL BOXER(II,JJ,XOUT,YOUT,DOVER) IF(DOVER.GT.0.0) THEN C Re-normalise the area overlap using the Jacobian DOVER=DOVER/JACO C Count the hits NHIT=NHIT+1 VC=NCOU(II,JJ) C Allow for stretching because of scale change D=DATA(I)*S2 C Scale the weighting mask by the scale factor C and inversely by the Jacobian to ensure conservation C of weight in the output W=WEI(I)*WTSCL IF(W.NE.1.0) THEN IF(VC.EQ.0.0) THEN NDAT(II,JJ)=D NCOU(II,JJ)=DOVER*W ELSE NDAT(II,JJ)= : (NDAT(II,JJ)*VC+DOVER*W*D)/ : (VC+DOVER*W) NCOU(II,JJ)=VC+DOVER*W ENDIF ELSE IF(VC.EQ.0.0) THEN NDAT(II,JJ)=D NCOU(II,JJ)=DOVER ELSE NDAT(II,JJ)= : (NDAT(II,JJ)*VC+DOVER*D)/ : (VC+DOVER) NCOU(II,JJ)=VC+DOVER ENDIF ENDIF ENDIF ENDIF ENDDO ENDDO C Count cases where the pixel is off the output image IF(NHIT.EQ.0) NMISS=NMISS+1 ENDDO ENDDO C Report number of points which were outside the output frame IF(NMISS.GT.0) THEN WRITE(CHARS,'(''! Warning, '',I7,'' points were '','// : '''outside the output image.'')') : NMISS CALL UMSPUT(CHARS,1,0,ISTAT) ENDIF C We now try to update the WCS C First CRPIX1,CRPIX2,CRVAL1,CRVAL2 IF(NONLIN) THEN XNL=EVALCU(REAL(WCS(1)-XCEN),REAL(WCS(3)-YCEN),XCO) YNL=EVALCU(REAL(WCS(1)-XCEN),REAL(WCS(3)-YCEN),YCO) XNLD=EVALCU(REAL(WCS(1)-XCEN)+1.0,REAL(WCS(3)-YCEN),XCO) YNLD=EVALCU(REAL(WCS(1)-XCEN),REAL(WCS(3)-YCEN)+1.0,YCO) ELSE XNL=WCS(1)-XCEN YNL=WCS(3)-YCEN XNLD=XNL+1 YNLD=YNL+1 ENDIF C Now the centre position IF(ROTFIR) THEN WCS(1)=DBLE(XC*XNL-YS*YNL+XT) WCS(3)=DBLE(XS*XNL+YC*YNL+YT) ELSE WCS(1)=DBLE(XC*(XNL+XSH)-YS*(YNL+YSH)+XP) WCS(3)=DBLE(XS*(XNL+XSH)+YC*(YNL+YSH)+YP) ENDIF C The reference point on the sky is unchanged (WCS(2),WCS(4)) C Now we scale and rotate the CD matrix (WCS(5-8)) C First calculate the determinant for the 2x2 matrix inversion DET=COSTH*(XNLD-XNL)*(YNLD-YNL)*COSTH+ : SINTH*(XNLD-XNL)*(YNLD-YNL)*SINTH A=COSTH*(YNLD-YNL)*WCS(5)-SINTH*(XNLD-XNL)*WCS(7) B=COSTH*(YNLD-YNL)*WCS(6)-SINTH*(XNLD-XNL)*WCS(8) C=SINTH*(YNLD-YNL)*WCS(5)+COSTH*(XNLD-XNL)*WCS(7) D=SINTH*(YNLD-YNL)*WCS(6)+COSTH*(XNLD-XNL)*WCS(8) WCS(5)=DBLE(A/XF/DET) WCS(6)=DBLE(B/XF/DET) WCS(7)=DBLE(C/XF/DET) WCS(8)=DBLE(D/XF/DET) C Set good status if we get this far ISTAT=0 RETURN END C C BOXER -- compute area of box overlap C C Calculate the area common to input clockwise polygon x(n), y(n) with C square (is, js) to (is+1, js+1). C This version is for a quadrilateral. C C W.B. Sparks STScI 2-June-1990. C Phil Hodge 20-Nov-1990 Change calling sequence; single precision. C Richard Hook ECF 24-Apr-1996 Change coordinate origin C so that centre of pixel has integer position SUBROUTINE BOXER (IS, JS, X, Y, DAREA) IMPLICIT NONE INTEGER IS, JS REAL X(*), Y(*) REAL DAREA C-- REAL PX(4), PY(4), SUM REAL SGAREA INTEGER I C Set up coords relative to unit square at origin C Note that the +0.5s were added when this code was C included in DRIZZLE DO I = 1, 4 PX(I) = X(I) - IS +0.5 PY(I) = Y(I) - JS +0.5 ENDDO C For each line in the polygon (or at this stage, input quadrilateral) C calculate the area common to the unit square (allow negative area for C subsequent `vector' addition of subareas). SUM = 0.0 DO I = 1, 3 SUM = SUM + SGAREA (PX(I), PY(I), PX(I+1), PY(I+1), IS, JS) ENDDO SUM = SUM + SGAREA (PX(4), PY(4), PX(1), PY(1), IS, JS) DAREA = SUM RETURN END REAL FUNCTION SGAREA (X1, Y1, X2, Y2, IS, JS) C C To calculate area under a line segment within unit square at origin. C This is used by BOXER C IMPLICIT NONE REAL X1, Y1, X2, Y2 INTEGER ISTAT,IS,JS REAL M, C, DX REAL XLO, XHI, YLO, YHI, XTOP CHARACTER*80 CHARS LOGICAL NEGDX DX = X2 - X1 C Trap vertical line IF (DX .EQ. 0.0) THEN SGAREA = 0.0 GO TO 80 ENDIF C Order the two input points in x IF (X1 .LT. X2) THEN XLO = X1 XHI = X2 ELSE XLO = X2 XHI = X1 ENDIF C And determine the bounds ignoring y for now IF (XLO .GE. 1.0) THEN SGAREA = 0.0 GO TO 80 ENDIF IF (XHI .LE. 0.0) THEN SGAREA = 0.0 GO TO 80 ENDIF XLO = MAX (XLO, 0.0) XHI = MIN (XHI, 1.0) C Now look at y C basic info about the line y = mx + c NEGDX = (DX .LT. 0.0) M = (Y2 - Y1) / DX C = Y1 - M * X1 YLO = M * XLO + C YHI = M * XHI + C C Trap segment entirely below axis IF (YLO .LE. 0.0 .AND. YHI .LE. 0.0) THEN SGAREA = 0.0 GO TO 80 ENDIF C Adjust bounds if segment crosses axis (to exclude anything below axis) IF (YLO .LT. 0.0) THEN YLO = 0.0 XLO = -C/M ENDIF IF (YHI .LT. 0.0) THEN YHI = 0.0 XHI = -C/M ENDIF C There are four possibilities: both y below 1, both y above 1 C and one of each. IF (YLO .GE. 1.0 .AND. YHI .GE. 1.0) THEN C Line segment is entirely above square IF (NEGDX) THEN SGAREA = XLO - XHI ELSE SGAREA = XHI - XLO ENDIF GO TO 80 ENDIF IF (YLO .LE. 1.0 .AND. YHI .LE. 1.0) THEN C Segment is entirely within square IF (NEGDX) THEN SGAREA = 0.5 * (XLO-XHI) * (YHI+YLO) ELSE SGAREA = 0.5 * (XHI-XLO) * (YHI+YLO) END IF GO TO 80 ENDIF C otherwise it must cross the top of the square XTOP = (1.0 - C) / M C Check for a problem - note that this now has a small C error margin included to avoid errors coming from numerical C inaccuracy IF (XTOP .LT. XLO-1.0E-5 .OR. XTOP .GT. XHI+1.0E-5) THEN WRITE(CHARS,'(''! Warning, box overlap calculation'','// : ''' problem at ('',I5,'','',I5,'')'')') IS,JS CALL UMSPUT(CHARS,1,0,ISTAT) ENDIF IF (YLO .LT. 1.0) THEN IF (NEGDX) THEN SGAREA = -(0.5 * (XTOP-XLO) * (1.0+YLO) + XHI - XTOP) ELSE SGAREA = 0.5 * (XTOP-XLO) * (1.0+YLO) + XHI - XTOP END IF GO TO 80 ENDIF IF (NEGDX) THEN SGAREA = -(0.5 * (XHI-XTOP) * (1.0+YHI) + XTOP-XLO) ELSE SGAREA = 0.5 * (XHI-XTOP) * (1.0+YHI) + XTOP-XLO END IF 80 CONTINUE RETURN END