C @(#)rebdec.for 17.1.1.1 (ESO-DMD) 01/25/02 17:19:21 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM REBDEC C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program REBDEC version 1.00 880930 C Basic algorithm: L. Lucy ESO - Garching C MIDAS implementation: D. Baade ESO - Garching C K. Banse 1.10 890517 C 1.20 900220 C C.KEYWORDS C (nonlinear) rebinning, deconvolution, point spread function, undersampling C C.PURPOSE C Rebin an image to smaller stepsize to remove effects of coarse PSF sampling C ("pixelation"). Optionally deconvolve with PSD. C C.ALGORITHM C Get the names of input + output frames from IN_A + OUT_A (both 2-D) C Do the rebinning (which is nonlinear in the way the flux is redistributed C whereas the spatial coordinate transformation is strictly linear) while C simultaneously deconvolving with user-supplied point spread function. C C.INPUT/OUTPUT C the following keys are used: C C IN_A/C/1/60 input frame C IN_B/C/1/60 point spread function C OUT_A/C/1/60 output frame C INPUTI/I/1/2 zoom factors in x and y C INPUTI/I/3/1 number of iterations C C.VERSIONS C 1.10 convert to Portable Midas C 1.20 fix bug in STIPUT (missing NAXIS..) C C 010201 last modif C C -------------------------------------------------------------------- C IMPLICIT NONE C INTEGER STAT,ZOOM(2),NITER,NAXIN,NPXIN(2) INTEGER INNO,PSFNO,OUTNO,DUMNO,D0NO,D1NO,D2NO,D3NO,D4NO INTEGER NAXPSF,NPXPSF(2),IAV,OFF(2) INTEGER NPXPSR(4),NAXPSR,I,SUBLO(2),GS INTEGER NPXOUT(2),RESPIX(2),SIZE,NAXOUT,NPIXL(2) INTEGER NPIXZ(2),MSIZ,INPUTI(4) INTEGER UNI(1),NULO,MADRID(1) INTEGER*8 INPTR,PSFPTR,OUTPTR,DUMPTR,D0PTR,D1PTR,D2PTR, + D3PTR,D4PTR C CHARACTER INFRM*60,OUTFRM*60,PSF*60 CHARACTER CUNIN*64,CUNPSF*64,IDENT*72 CHARACTER IDENTP*72,CUNOUT*64,IDENTO*72 C DOUBLE PRECISION STAIN(2),STPIN(2),STAPSF(2) DOUBLE PRECISION STPPSF(2),STPOUT(2),STPPSR(4) DOUBLE PRECISION STAPSR(4),STAOUT(2),TEST,EPS C REAL CUTS(2),CUTVALS(4) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C COMMON /VMR/ MADRID C DATA CUTVALS /4*0./, SUBLO /1,1/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C set up MIDAS environment + enable automatic error abort CALL STSPRO('REBDEC') C C get name of input image + map frame CALL STKRDC('IN_A',1,1,60,IAV,INFRM,UNI,NULO,STAT) CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXIN,NPXIN,STAIN,STPIN, + IDENT,CUNIN,INPTR,INNO,STAT) IF (NAXIN.NE.2) CALL STETER + (1,'ERROR: Currently only 2-D frames supported!') C C get and check zoom factors in x + y ... and the number of iterations CALL STKRDI('INPUTI',1,4,IAV,INPUTI,UNI,NULO,STAT) ZOOM(1) = INPUTI(1) ZOOM(2) = INPUTI(2) NITER = INPUTI(3) GS = INPUTI(4) C IF ((((ZOOM(1)*ZOOM(2))/2)*2).EQ.(ZOOM(1)*ZOOM(2))) + CALL STETER + (2,'ERROR: Both zoom factors need to be odd!') IF ((ZOOM(1).LT.1).OR.(ZOOM(2).LT.1)) + CALL STETER + (3,'ERROR: Zoom factors less than unity are illegal!') C IF (NITER.LT.1) CALL STETER + (4,'ERROR: Number of iterations must at least be 1!') C C get name of point spread function frame + map it CALL STKRDC('IN_B',1,1,60,IAV,PSF,UNI,NULO,STAT) CALL STIGET(PSF,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXPSF,NPXPSF,STAPSF, + STPPSF,IDENTP,CUNPSF,PSFPTR,PSFNO,STAT) C C Further error traps: make sure dimensions are the same ... IF (NAXPSF.NE.NAXIN) + CALL STETER + (5,'Dimensions of inframe and PSF must be the same!') C C ... steps have the right relative proportions, IF ((NINT(STPIN(1)/STPPSF(1)).NE.ZOOM(1)).OR. + (NINT(STPIN(2)/STPPSF(2)).NE.ZOOM(2))) + CALL STETER + (6,'Steps of inframe and PSF must differ by zoom fact!') C C ... the zoom factor is not even when the number of pixels in the input PSF C is odd (check both X and Y) DO 100 I=1,2 IF ((NPXPSF(I)/2)*2.EQ.NPXPSF(I)) + CALL STETER + (7,'NPIX(PSF) must be odd in both X and Y!') 100 CONTINUE C C Get name of output frame, ... CALL STKRDC('OUT_A',1,1,60,IAV,OUTFRM,UNI,NULO,STAT) IF (OUTFRM.EQ.INFRM) ! in- and output frame must be different + CALL STETER + (8,'Names of inframe + outframe must be different!') C C ... depending on its existence or not, C C ... open it for reading and writing or create it for writing to it C IF (GS.EQ.0) THEN ! we start without an initial guess C C Calculate further defining parameters of output frame ... STPOUT(1)=STPIN(1)/ZOOM(1) STPOUT(2)=STPIN(2)/ZOOM(2) NPXOUT(1)=NPXIN(1)*ZOOM(1) NPXOUT(2)=NPXIN(2)*ZOOM(2) C C The zooming STAOUT(1)=STAIN(1)-STPOUT(1)*(ZOOM(1)/2) STAOUT(2)=STAIN(2)-STPOUT(2)*(ZOOM(2)/2) C C ... and map it: IDENTO(1:) = IDENT(1:) CUNOUT(1:) = CUNIN(1:) NAXOUT = 2 CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXOUT,NPXOUT,STAOUT, + STPOUT,IDENT,CUNIN,OUTPTR,OUTNO,STAT) C ELSE ! there is a guess C C If there is a guess, we first map the frame ... CALL STIGET(OUTFRM,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE, + 2,NAXOUT,NPXOUT,STAOUT, + STPOUT,IDENTO,CUNOUT,OUTPTR,OUTNO,STAT) C C ... and then check if it is technically compatible with the C input frame and the zoom factors: C IF (NAXOUT.NE.2) CALL STETER + (9,'Output frame, too, must be 2-dimensional!') IF ((NPXIN(1)*NPXIN(2)*ZOOM(1)*ZOOM(2)).NE. + (NPXOUT(1)*NPXOUT(2))) + CALL STETER + (10,'NPIX(OUT) incompatible with NPIX(IN) + zoom factors!') TEST = (STAOUT(1)+STPOUT(1)*(ZOOM(1)/2))* + (STAOUT(2)+STPOUT(2)*(ZOOM(2)/2)) EPS = STPIN(1)*STPIN(2)*0.01 IF (ABS(TEST-(STAIN(1)*STAIN(2))).GT.EPS) + CALL STETER + (11,'START values of in- and outframe differ!') EPS = STPIN(1)*STPIN(2)*0.0001 IF (ABS((STPOUT(1)*ZOOM(1)*STPOUT(2)*ZOOM(2))- + (STPIN(1)*STPIN(2))).GT.EPS) + CALL STETER + (12,'steps of in- and outframe must differ by factor ZOOM!') ENDIF ! guess or no guess? C C We also need six working buffers for ... C C a) the output frame with extrapolation beyond its edges DO 400 I=1,2 C C The radius of the PSF determines the width of the margin to be added C around inframe, for outframe additionally the zoom factor must be considered: OFF(I) = NPXPSF(I)/(2*ZOOM(I)) + 2 NPIXL(I) = NPXIN(I)+2*OFF(I) NPIXZ(I) = NPIXL(I)*ZOOM(I) 400 CONTINUE C MSIZ = NPIXZ(1) * NPIXZ(2) CALL STFCRE('DUMMY',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + MSIZ,DUMNO,STAT) CALL STFMAP(DUMNO,F_X_MODE,1,MSIZ,IAV,DUMPTR,STAT) C C If there is a guess, we transfer it right away to the enlarged frame: IF (GS.EQ.1) + CALL EXTND(MADRID(OUTPTR),MADRID(DUMPTR), + NPXOUT(1),NPXOUT(2),NPIXZ(1),NPIXZ(2),OFF, + ZOOM(1),ZOOM(2)) C C In order not to waste virtual memory space, we temporarily close OUTFRM CALL STFCLO(OUTNO,STAT) C C b) the input frame after extrapolation beyond its edges MSIZ = NPIXL(1) * NPIXL(2) CALL STFCRE('DUMMY0',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + MSIZ,D0NO,STAT) CALL STFMAP(D0NO,F_X_MODE,1,MSIZ,IAV,D0PTR,STAT) C C We do the enlargement here and then close inframe (no longer needed) C CALL EXTND (MADRID(INPTR),MADRID(D0PTR),NPXIN(1),NPXIN(2), + NPIXL(1),NPIXL(2),OFF,1,1) CALL STFCLO(INNO,STAT) C C c) ratio of observed to calculated image MSIZ = NPIXL(1) * NPIXL(2) CALL STFCRE('DUMMY1',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + MSIZ,D1NO,STAT) CALL STFMAP(D1NO,F_X_MODE,1,MSIZ,IAV,D1PTR,STAT) C C d) correction factors to be applied to previous intermediate result MSIZ = NPIXZ(1) * NPIXZ(2) CALL STFCRE('DUMMY2',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + MSIZ,D2NO,STAT) CALL STFMAP(D2NO,F_X_MODE,1,MSIZ,IAV,D2PTR,STAT) C C e) re-calculated frame on coarse grid MSIZ = NPIXL(1) * NPIXL(2) CALL STFCRE('DUMMY3',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + MSIZ,D3NO,STAT) CALL STFMAP(D3NO,F_X_MODE,1,MSIZ,IAV,D3PTR,STAT) C C f) PSF, but re-sampled on coarse grid for various centerings. C C First, the dimension of this new PSF need to be properly determined: NAXPSR = 4 DO 1000 I=1,2 C C The number of pixels of the coarse PSF (see subroutine PSFRESAMP) must be C odd. If ZOOM divides NPIX(PSF), this is automatically fulfilled because C both numers are odd. NPXPSR(2*I-1) = NPXPSF(I) / ZOOM(I) C C Otherwise, we just make it odd: IF (((NPXPSF(I)/ZOOM(I))*ZOOM(I)).NE.NPXPSF(I)) + NPXPSR(2*I-1) = (NPXPSR(2*I-1)/2)*2 + 1 C C Finally, an additional margin of 1 is added at either limit in order C to avoid truncations of input-PSF: NPXPSR(2*I-1) = NPXPSR(2*I-1)+2 NPXPSR(2*I) = ZOOM(I) STPPSR(2*I-1) = STPPSF(I) STPPSR(2*I) = 1.D0/ZOOM(I) STAPSR(2*I-1) = STAPSF(I) STAPSR(2*I) = 1.D0/ZOOM(I) 1000 CONTINUE C MSIZ = NPXPSR(1) * NPXPSR(2) * NPXPSR(3) * NPXPSR(4) CALL STFCRE('DUMMY4',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, + MSIZ,D4NO,STAT) CALL STFMAP(D4NO,F_X_MODE,1,MSIZ,IAV,D4PTR,STAT) C C Do the resampling of the PSF here and then close the input PSF C because it's no longer used CALL PSFRES(NPXPSF(1),NPXPSF(2),NPXPSR(1),NPXPSR(3), + ZOOM(1),ZOOM(2),MADRID(PSFPTR),MADRID(D4PTR)) CALL STFCLO(PSFNO,STAT) C C Now comes the real thing: CALL REBDC(NPIXL(1),NPIXL(2),NPIXZ(1),NPIXZ(2),NITER, + NPXPSR(1),NPXPSR(3),ZOOM(1),ZOOM(2),GS, + MADRID(D0PTR),MADRID(D1PTR),MADRID(D2PTR), + MADRID(D3PTR),MADRID(D4PTR),MADRID(DUMPTR)) C C Close dispensible intermediate files and re-open outframe . Cut off C extrapolated regions of intermediate result as to match size of OUTFRM, C insert into outframe: C CALL STFCLO(D0NO,STAT) CALL STFCLO(D1NO,STAT) CALL STFCLO(D2NO,STAT) CALL STFCLO(D3NO,STAT) CALL STFCLO(D4NO,STAT) C CALL STIGET(OUTFRM,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE, + 2,NAXOUT,NPXOUT,STAOUT, + STPOUT,IDENTO,CUNOUT,OUTPTR,OUTNO,STAT) CALL TRUNC(MADRID(DUMPTR),MADRID(OUTPTR),NPIXZ(1),NPIXZ(2), + NPXOUT(1),NPXOUT(2),OFF,ZOOM(1),ZOOM(2)) C C Last, we determine physical minimum and maximum of OUTFRM CALL STVALS('MIN',MADRID(OUTPTR),2,NPXOUT,SUBLO, + NPXOUT,CUTVALS,CUTS,RESPIX,SIZE,STAT) C C Copy descriptors and update LHCUTS + HISTORY CALL DSCUPT(INNO,OUTNO,' ',STAT) CALL STDWRR(OUTNO,'LHCUTS',CUTS,3,2,UNI,STAT) C C Bye, Bye CALL STSEPI END SUBROUTINE REBDC(ICX,ICY,I1X,I1Y,NITER,OPSFX,OPSFY,ZX,ZY,GS, 1 PHT,RT,CF,PHI,PSF,PSI) C C Removal of pixel grid from CCD data plus seeing by appropriate C rebinning (including 2-D deconvolution) to smaller step size. C The deconvolution algorithm follows the description given in C The Astronomical Journal, Vol. 79, p. 745-754 (1974) but the C rebinning to a smaller step is an extra feature. C C L.B. Lucy ESO - Garching 1988 August C Modifications applied for MIDAS version (incl. generalisation C to non-square images and PSF's): C C D. Baade ESO - Garching 880930 C IMPLICIT NONE C INTEGER ICX,ICY,GS,ZX,ZY,NITER,STAT,I1X,I1Y,IT,I,J INTEGER OPSFX,OPSFY C REAL PHT(ICX,ICY),PHI(ICX,ICY),CD,CD1,RT(ICX,ICY) REAL PSF(OPSFX,ZX,OPSFY,ZY),CF(I1X,I1Y),PSI(I1X,I1Y) C CHARACTER TIME*25,CBUF*80 C C PHT inframe, but extrapolated beyond edges, dims: ICX x ICY C PHI calculated frame in coarse grid C RT basically, the ratio observed/calculated in coarse grid C CF correction factor to be applied to current fine-grid result C PSI fine-grid output frame, STAINg from some trivial guess C IC dimensions of input frame: ICX x ICY C Z zoom factors coarse -> fine grid: ZX, ZY C I1 dimensions of output frame: ZX*ICX x ZY*ICY C C C ... If an initial guess has not been provided, we start with a C ... smooth trial-image C CD = ZX*ZY IF (GS.EQ.0) THEN CD1 = CD/(I1X*I1Y) DO 200 I=1,I1X DO 100 J=1,I1Y PSI(I,J) = CD1 ! initial smooth trial - image 100 CONTINUE 200 CONTINUE ENDIF ! start frame prepared C DO 1000 IT=1,NITER ! iteration loop C C Form PSF-mapped coarse-grid image of guess PSI (result in frame PHI): CALL INTEGD(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY, + PSF,PSI,PHI) C C Compute ratio of observed to calculated (guessed) image C DO 400 I=1,ICX DO 300 J=1,ICY RT(I,J) = 1.0 ! initialize IF (PHI(I,J).NE.0.0) + RT(I,J) = CD*PHT(I,J)/PHI(I,J) ! obs/calc 300 CONTINUE 400 CONTINUE C C Determine correction factor to current intermediate result CALL INTEGF(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY, + PSF,RT,CF) C C Apply the correction, thereby completing this iteration DO 600 I=1,I1X DO 500 J=1,I1Y PSI(I,J) = PSI(I,J)*CF(I,J) ! apply correc factor 500 CONTINUE 600 CONTINUE C C The code is slow. Progress messages seem a good idea: CALL GENTIM(TIME) WRITE (CBUF,10000) TIME,IT CALL STTPUT(CBUF,STAT) C 1000 CONTINUE C RETURN C 10000 FORMAT (A,' ',I2,'. iteration completed') END SUBROUTINE EXTND (RAW,PHT,IRX,IRY,ICX,ICY,OFF,ZX,ZY) C C To avoid ugly effects at the frame limits, repeat the outermost C columns (rows) RX (RY) times, i.e. by the halfwidth of the C PSF resampled on the coarse grid. C IMPLICIT NONE C INTEGER IRX,IRY,ICX,ICY,I,J,IQ,JQ,OFF(2),ZX,ZY,OFX,OFY C REAL RAW(IRX,IRY),PHT(ICX,ICY) C OFX = OFF(1)*ZX OFY = OFF(2)*ZY C C Insert raw frame at center of output frame C DO 200 I=1,IRX IQ = I+OFX DO 100 J=1,IRY JQ = J+OFY PHT(IQ,JQ) = RAW(I,J) 100 CONTINUE 200 CONTINUE C C ... Do the extrapolation, first in X: C DO 500 I=1,OFX DO 300 J=1,IRY PHT(I,J+OFY) = RAW(1,J) PHT(IRX+OFX+I,J+OFY) = RAW(IRX,J) 300 CONTINUE C C Here we also deal with the four corners C DO 400 J=1,OFY PHT(I,J) = RAW(1,1) PHT(IRX+OFX+I,J) = RAW(IRX,1) PHT(I,IRY+OFY+J) = RAW(1,IRY) PHT(IRX+OFX+I,IRY+OFY+J) = RAW(IRX,IRY) 400 CONTINUE 500 CONTINUE C C The continuation in Y: C DO 700 I=1,IRX DO 600 J=1,OFY PHT(I+OFX,J) = RAW(I,1) PHT(I+OFX,IRY+OFY+J) = RAW(I,IRY) 600 CONTINUE 700 CONTINUE C RETURN END SUBROUTINE PSFRES (IPSFX,IPSFY,OPSFX,OPSFY,ZX,ZY,PSFIN,PSF) C C Author: L. Lucy ESO - Garching 1988 August C Modifications: D. Baade ESO - Garching 880930 C C The subroutine resamples the user-supplied high-resolution PSF C on the coarser grid of the image to be deconvolved. The array C holding the resulting PSF is four-dimensional because the resampling C is done ZX x ZY times where each time the input PSF is C shifted by 1/ZX and/or 1/ZY as to cover all cases of the C relative positioning of the two PSF grids. C IMPLICIT NONE C INTEGER II,JJ,IP,JP,K,L,ZX,ZY,INDX,INDY,IPSFX,IPSFY INTEGER OPSFX,OPSFY,DX,DY,OFFX,OFFY,CX,CY C REAL PSFIN(IPSFX,IPSFY),PSFBUF,CD,PSF(OPSFX,ZX,OPSFY,ZY) C CD = 1./(ZX*ZY) DX = (OPSFX*ZX-IPSFX)/2 ! the excess in width of the resampled PSF DY = (OPSFY*ZY-IPSFY)/2 ! is always an even number of small pixels. C C Shift by one small pixel in X ... C DO 1000 II=1,ZX OFFX = II-(ZX+1)/2+DX C C and/or Y DO 800 JJ=1,ZY OFFY = JJ-(ZY+1)/2+DY C C For all coarse-grid pixels ... DO 600 IP=1,OPSFX CX = (IP-1)*ZX DO 500 JP=1,OPSFY CY = (JP-1)*ZY PSFBUF = 0. ! initialize C C co-add all fine-grid pixels covered by them at the current offsets DO 400 K=1,ZX INDX = CX+K-OFFX C C Skip fine-grid pixels for which the input-PSF is undefined: IF ( (INDX.LT.1) .OR. (INDX.GT.IPSFX) ) GOTO 400 DO 300 L=1,ZY INDY = CY+L-OFFY IF ( (INDY.GE.1) .AND. (INDY.LE.IPSFY) ) + PSFBUF = PSFBUF + PSFIN(INDX,INDY) 300 CONTINUE 400 CONTINUE C PSF(IP,II,JP,JJ) = PSFBUF*CD 500 CONTINUE 600 CONTINUE C 800 CONTINUE 1000 CONTINUE C RETURN END SUBROUTINE INTEGD(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY, + PSF,PSI,PHI) C C Author: L. Lucy ESO - Garching 1988 August C Modifications: D. Baade ESO - Garching 880930 C C This subroutine maps the current fine-grid guess, PSI, C on to the coarse grid (result in array PHI). The mapping C uses the stack of PSF's that have been constructed by C re-sampling the input PSF onto the coarse grid, i.e. basically C a convolution with the psf is done. C IMPLICIT NONE C INTEGER ICX,ICY,I1X,I1Y,OPSFX,OPSFY,IMIN,IMAX,JMIN,JMAX INTEGER ID,ID0,ID1,ID2,II,IP,JD,JD0,JD1,JD2,JP,JJ,ZX,ZY INTEGER I,J,CX(10000),CY(10000),RX,RY,RX1,RY1,RX1I,RY1J INTEGER AUXX(2000),AUXY(2000) C REAL PSF(OPSFX,ZX,OPSFY,ZY),PSI(I1X,I1Y),PHIBUF,PHI(ICX,ICY) C RX = OPSFX/2 ! 'radius' of re-sampled PSF, RY = OPSFY/2 ! full diameter is 2*R+1 RX1 = RX+1 RY1 = RY+1 C C For all fine pixels, determine X and Y coordinates on the coarse grid. C Also, introduce auxiliary array to further save computing time. C DO 100 I=1,I1X CX(I) = (I-1)/ZX + 1 AUXX(I) = (CX(I)-1)*ZX 100 CONTINUE DO 200 J=1,I1Y CY(J) = (J-1)/ZY + 1 AUXY(J) = (CY(J)-1)*ZY 200 CONTINUE C C Proceed pixel by pixel on the coarse grid C IMIN = RX1 IMAX = ICX-RX JMIN = RY1 JMAX = ICY-RY ID0 = 1-RX1*ZX JD0 = 1-RY1*ZY C DO 2000 I=IMIN,IMAX ! pixel of coarse grid ID1=ID0+I*ZX ID2=(I+RX)*ZX RX1I=RX1-I DO 1800 J=JMIN,JMAX ! pixel of coarse grid JD1=JD0+J*ZY JD2=(J+RY)*ZY RY1J=RY1-J C C Then sum over all relevant pixels in fine grid. the PSF serves as C the weighting function. For a small pixel (II,JJ) within a given large C pixel (I,J) the value PSF(IP,II,JP,JJ) is taken. IP and JP point C to the coarse pixels contributing (via the PSF) to the currently C treated coarse pixel. IP and JP are no absolute coordinates but C rather are the distances of the contributing large pixels from the C large pixel whose flux is just being computed. II and JJ are also C relative indices and only range from 1 to ZX and ZY, respectively. C PHIBUF = 0. ! initialize C DO 1600 ID=ID1,ID2 ! relevant range in fine grid II=ID-AUXX(ID) ! position within large pixel IP=CX(ID)+RX1I ! position of large pixel in PSF DO 1500 JD=JD1,JD2 ! relevant range in fine grid JJ = JD-AUXY(JD) ! position within large pixel JP = CY(JD)+RY1J ! position of large pixel in PSF PHIBUF = PHIBUF+PSF(IP,II,JP,JJ)*PSI(ID,JD) ! off we go 1500 CONTINUE 1600 CONTINUE PHI(I,J) = PHIBUF ! this coarse bin finished C 1800 CONTINUE 2000 CONTINUE C C Now extrapolate solution achieved so far to edges of frame by C constant extrapolation C DO 3000 I=1,RX DO 2200 J=JMIN,JMAX PHI(I,J)=PHI(IMIN,J) PHI(IMAX+I,J)=PHI(IMAX,J) 2200 CONTINUE DO 2400 J=1,RY PHI(I,J)=PHI(IMIN,JMIN) PHI(IMAX+I,J)=PHI(IMAX,JMIN) PHI(I,JMAX+J)=PHI(IMIN,JMAX) PHI(IMAX+I,JMAX+J)=PHI(IMAX,JMAX) 2400 CONTINUE 3000 CONTINUE DO 4000 I =IMIN,IMAX DO 3300 J=1,RY PHI(I,J)=PHI(I,JMIN) PHI(I,JMAX+J)=PHI(I,JMAX) 3300 CONTINUE 4000 CONTINUE C RETURN END SUBROUTINE INTEGF(I1X,I1Y,ICX,ICY,ZX,ZY,OPSFX,OPSFY,PSF,RT,CF) C C Author: L. Lucy ESO - Garching 1988 August C Modifications: D. Baade ESO - Garching 880930 C C Here the correction factors (pixelwise on the fine grid, array CF) C are computed as the PSF-mapped ratio (pixelwise on the coarse grid, C array RT) of the observed image to the result of the last iteration C after its re-sampling on the coarse grid. C IMPLICIT NONE C INTEGER I,J,I1X,I1Y,ICX,ICY,II,ID1,ID2,IMIN,IMAX,JMIN,JMAX INTEGER JJ,JD1,JD2,ID,IP,JD,JP,OPSFX,OPSFY,ZX,ZY,RX,RY INTEGER CX(10000),CY(10000),AUXX1(10000),AUXY1(10000) INTEGER AUXX2(10000),AUXY2(10000) C REAL PSF(OPSFX,ZX,OPSFY,ZY),CF(I1X,I1Y),RT(ICX,ICY) REAL CFBUF C RX=OPSFX/2 ! 'radius' of coarse PSF, RY=OPSFY/2 ! full diameter is 2*R+1 C IMIN=RX*ZX+1 IMAX=I1X-RX*ZX JMIN=RY*ZY+1 JMAX=I1Y-RY*ZY C C Determine for all fine pixels coordinates on coarse grid (taking C into account that the input frame has been enlarged). C Also, introduce auxiliary arrays to further reduce computing time. C DO 100 I=IMIN,IMAX CX(I)=(I-1)/ZX+1 ! X-coordinates on coarse grid AUXX1(I)=RX-CX(I)+1 AUXX2(I)=(CX(I)-1)*ZX 100 CONTINUE DO 200 J=JMIN,JMAX CY(J)=(J-1)/ZY+1 ! Y-coordinates on coarse grid AUXY1(J)=RY-CY(J)+1 AUXY2(J)=(CY(J)-1)*ZY 200 CONTINUE C DO 1000 I=IMIN,IMAX ! pixel on fine grid II=I-AUXX2(I) ! pixel-# within coarse pixel ID1=CX(I)-RX ID2=CX(I)+RX C DO 900 J=JMIN,JMAX ! pixel on fine grid JJ=J-AUXY2(J) ! pixel-# within coarse pixel JD1=CY(J)-RY JD2=CY(J)+RY C CFBUF = 0. ! initialize C DO 800 ID=ID1,ID2 ! relevant coarse-grid pixels IP=AUXX1(I)+ID ! X-position within coarse PSF DO 700 JD=JD1,JD2 ! relevant coarse-grid pixels JP=AUXY1(J)+JD ! Y-position within coarse PSF CFBUF=CFBUF+PSF(IP,II,JP,JJ)*RT(ID,JD) ! the correction 700 CONTINUE 800 CONTINUE CF(I,J) = CFBUF C C At this point the corrections for fine-grid pixel (I,J) to the C ratios at positions (ID,JD) on the coarse grid have been fully C evaluated for all pixels (ID,JD) concerned by pixel (I,J). C 900 CONTINUE 1000 CONTINUE C C Fill up rest of frame by constant extrapolation of the edges C DO 2000 I=1,RX*ZX DO 1100 J=JMIN,JMAX CF(I,J)=CF(IMIN,J) CF(IMAX+I,J)=CF(IMAX,J) 1100 CONTINUE DO 1200 J=1,RY*ZY CF(I,J)=CF(IMIN,JMIN) CF(IMAX+I,J)=CF(IMAX,JMIN) CF(I,JMAX+J)=CF(IMIN,JMAX) CF(IMAX+I,JMAX+J)=CF(IMAX,JMAX) 1200 CONTINUE 2000 CONTINUE DO 3000 I =IMIN,IMAX DO 2500 J=1,RY*ZY CF(I,J)=CF(I,JMIN) CF(I,JMAX+J)=CF(I,JMAX) 2500 CONTINUE 3000 CONTINUE C RETURN END SUBROUTINE TRUNC(PSI,OUT,I1X,I1Y,IOX,IOY,OFF,ZX,ZY) C C This subroutine truncates the result frame such as to correspond C to the true dimensions (in world coordinates) of the input frame. C IMPLICIT NONE C INTEGER I1X,I1Y,IOX,IOY,I,J,IQ,JQ,OFF(2),ZX,ZY,OFX,OFY C REAL PSI(I1X,I1Y),OUT(IOX,IOY) C OFX = OFF(1)*ZX OFY = OFF(2)*ZY C DO 100 I=1,IOX IQ=I+OFX DO 50 J=1,IOY JQ = J+OFY OUT(I,J) = PSI(IQ,JQ) 50 CONTINUE 100 CONTINUE C RETURN END