C @(#)deconv.for 13.1.1.1 (ES0-DMD) 06/02/98 18:17:32 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 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding 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 DECONV C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C program DECONV version 1.0 881028 C P. GROSBOL ESO - GARCHING C K. Banse (streamlined for MIDAS) C M Peron 900118, change the dimension of NPSF and the test for the C PSF centered at the origin C.KEYWORDS: C deconvolution C C.PURPOSE: C deconvolve an image with a point spread function C images may have up to 3 dimensions C C.ALGORITHM: C an iterative deconvolution method as published by L.B. Lucy C REF.: LUCY,.L.B 1974, A.J. VOL.79 , P.745 C C.INPUT/OUTPUT: C C the following keywords are used: 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 (1) = no. of iterations C if no_iter = 0, a convolution will be done ! C (2) = continue_flag C.VERSIONS C 1.00 from version 3.80 as of 870327 C use FORTRAN 77 + new ST interfaces C C------------------------------------------------------------ C IMPLICIT NONE C INTEGER I,IAV INTEGER*8 PNTRB,PNTRI,PNTRO,PNTRP,PSFPTR,PNTRX INTEGER IMNOI,IMNOO,IMNOP,IMNPSF INTEGER ISIZE,ISTAT INTEGER N,NAXI,NAXO,NAXP,NOI INTEGER NPIXI(3),NPIXO(3),NPIXP(3),INPUTI(2) INTEGER NPSF(3),LOW(3),OUTLOW(3),HIGH(3) INTEGER N0X,N0Y,MAXPSF INTEGER UNI(1),NULO,MADRID(1) C CHARACTER*60 INFRM,OUTFRM,PSFFRM CHARACTER CUNITI*64,CUNITO*64,CUNITP*64,MYACT*80 CHARACTER TEXT*72,IDENT*72 CHARACTER*80 ERMES3,ERMES4,ERMES5,ERMES6 CHARACTER TIME*40,MIDUNI*2,TMPFIL*14 C DOUBLE PRECISION STRI(3),STRO(3),STRP(3),T DOUBLE PRECISION STEPI(3),STEPO(3),STEPP(3),END C REAL CUTS(4) REAL RELC,EPSI REAL DECNVL !function C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/ MADRID C DATA CUTS /4*0./ DATA ERMES3 + /'Point Spread Function higher dimensioned than image '/ DATA ERMES4 + /'Point Spread Function not centered at origin (in x) '/ DATA ERMES6 + /'Point Spread Function not centered at origin (in y) '/ DATA ERMES5 /'Point Spread Function has different stepsize '/ DATA NPIXI /1,1,1/, NPIXO /1,1,1/, NPIXP /1,1,1/ DATA MAXPSF /10201/ !that's a 101*101 PSF DATA LOW /1,1,1/, OUTLOW /1,1,1/, HIGH /1,1,1/ DATA IDENT /' '/, CUNITI /' '/, CUNITP /' '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C get into MIDAS environment CALL STSPRO('DECONV') CALL STKRDC('MID$SESS',1,11,2,IAV,MIDUNI,UNI,NULO,ISTAT) TMPFIL(1:) = 'deco'//MIDUNI//'dum.bdf ' C C get input image and check dimensions CALL STKRDC('IN_A',1,1,60,IAV,INFRM,UNI,NULO,ISTAT) CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXI,NPIXI,STRI,STEPI,IDENT,CUNITI, + PNTRI,IMNOI,ISTAT) IF (NAXI.GT.2) NAXI = 2 !force to max. 2-dim frame C C get PSF and check it CALL STKRDC('IN_B',1,1,60,IAV,PSFFRM,UNI,NULO,ISTAT) CALL STIGET(PSFFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXP,NPIXP,STRP,STEPP,TEXT,CUNITP, + PNTRP,IMNOP,ISTAT) IF (NAXP.GT.NAXI) CALL STETER(1,ERMES3) C C handle x-dimension IF ( ABS(STEPP(1)-STEPI(1)) .GT. 0.1E-5*ABS(STEPI(1)) ) + CALL STETER(2,ERMES5) C END = STRP(1) + ( (NPIXP(1) - 1) * STEPP(1) ) T = ABS(STRP(1))-ABS(END) EPSI = 10.E-7 * ABS(STEPP(1)) IF (ABS(T).GT.EPSI)CALL STETER(3,ERMES4) T = MIN(ABS(STRP(1)),ABS(END)) C NPSF(1) = INT(T / ABS(STEPP(1))) !no. of pixels to left of center N0X = NINT( - (STRP(1)/STEPP(1)) ) + 1 !pixel no. of 0.0 C C for 1-dim files we are done already IF ( (NAXI.EQ.1).OR.(NAXP.EQ.1) ) THEN NPSF(2) = 1 LOW(1) = N0X - NPSF(1) !starting pixel in x HIGH(1) = N0X + NPSF(1) GOTO 2000 ENDIF C C handle y-dimension IF ( ABS(STEPP(2)-STEPI(2)) .GT. 0.1E-5*ABS(STEPI(2)) ) + CALL STETER(4,ERMES5) C END = STRP(2) + ( (NPIXP(2) - 1) * STEPP(2) ) T = ABS(STRP(2))-ABS(END) EPSI = 10.E-7 * ABS(STEPP(2)) IF (ABS(T).GT.EPSI)CALL STETER(5,ERMES6) T = MIN(ABS(STRP(2)),ABS(END)) C NPSF(2) = INT(T / ABS(STEPP(2))) !no. of pixels to left of center N0Y = NINT( - (STRP(2)/STEPP(2)) ) + 1 !pixel no. of 0.0 C C make minimum square PSF NPSF(1) = MIN(NPSF(1),NPSF(2)) NPSF(2) = NPSF(1) LOW(1) = N0X - NPSF(1) !starting pixel in x HIGH(1) = N0X + NPSF(1) LOW(2) = N0Y - NPSF(2) !starting pixel in y HIGH(2) = N0Y + NPSF(2) C C allocate virtual memory for working area = one line 2000 CALL STFXMP(NPIXI(1),D_R4_FORMAT,PNTRX,ISTAT) C C create PSF temporary STRP(1) = STRP(1) + ( (LOW(1)-1)*STEPP(1) ) STRP(2) = STRP(2) + ( (LOW(2)-1)*STEPP(2) ) NPSF(1) = HIGH(1) - LOW(1) + 1 !now really the dimensions... NPSF(2) = HIGH(2) - LOW(2) + 1 NPSF(3) = 1 MYACT(1:) = 'centered + normalized PSF ' CALL STIPUT(TMPFIL,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXP,NPSF,STRP,STEPP,MYACT, + CUNITP,PSFPTR,IMNPSF,ISTAT) CALL COPWND(MADRID(PNTRP),NPIXP,MADRID(PSFPTR),NPSF, + LOW,OUTLOW,HIGH) ISIZE = NPSF(1) * NPSF(2) IF (ISIZE.GT.MAXPSF) + CALL STTPUT('more than 10201 pixels in PSF...!',ISTAT) C C normalize termporary PSF CALL NORMFR(MADRID(PSFPTR),MADRID(PSFPTR),ISIZE) C C get result frame + check for continued iteration CALL STKRDC('OUT_A',1,1,60,IAV,OUTFRM,UNI,NULO,ISTAT) CALL STKRDI('INPUTI',1,2,IAV,INPUTI,UNI,NULO,ISTAT) NOI = INPUTI(1) IF (INPUTI(2).NE.1) THEN !iter_cont flag set ? CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXI,NPIXI,STRI,STEPI,IDENT,CUNITI, + PNTRO,IMNOO,ISTAT) CALL STDWRR(IMNOO,'LHCUTS',CUTS,1,4,UNI,ISTAT) CALL INITFR(MADRID(PNTRI),NAXI,NPIXI, + MADRID(PSFPTR),NAXP,NPSF, + MADRID(PNTRO),MADRID(PNTRX)) ELSE CALL STIGET(OUTFRM,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE, + 3,NAXO,NPIXO,STRO,STEPO,IDENT, + CUNITO,PNTRO,IMNOO,ISTAT) ENDIF C C update descriptors of result frame CALL DSCUPT(IMNOI,IMNOO,' ',ISTAT) IF (NOI.LT.0) GOTO 9000 !leave after one convolution C C create a temporary buffer frame ISIZE = 1 DO 5050, I=1,NAXI ISIZE = ISIZE * NPIXI(I) 5050 CONTINUE CALL STFXMP(ISIZE,D_R4_FORMAT,PNTRB,ISTAT) C C do the deconvolution + show processing times DO 6000, I=1,NOI RELC = DECNVL(MADRID(PNTRI),NAXI,NPIXI, + MADRID(PSFPTR),NAXP,NPIXP, + MADRID(PNTRO),MADRID(PNTRB), + MADRID(PNTRX)) CALL GENTIM(TIME) DO 5100, N=40,2,-1 IF (TIME(N:N).NE.' ') THEN IAV = N GOTO 5200 ENDIF 5100 CONTINUE IAV = 1 5200 WRITE(TEXT,10000) I,RELC MYACT(1:) = TEXT(1:36)//TIME(1:IAV)//' ' CALL STTPUT(MYACT,ISTAT) 6000 CONTINUE C C That's it folks... C 9000 CALL STSEPI C 10000 FORMAT('no.',I3,', rel. change: ',F10.3) END SUBROUTINE NORMFR(F,G,NO) C IMPLICIT NONE C INTEGER NO INTEGER I C REAL F(*),G(*) REAL R C DOUBLE PRECISION SUM C SUM = 0.D0 C DO 200 I=1,NO SUM = SUM + F(I) 200 CONTINUE C IF (ABS(SUM) .GT. 10.D-30) THEN R = SNGL( 1.D0/SUM ) DO 400 I=1,NO G(I) = F(I) * R 400 CONTINUE ELSE CALL STTPUT('Warning: Total flux = 0.0 !!!',I) ENDIF C RETURN END SUBROUTINE INITFR(F,NAF,NPF,PSF,NAP,NPP,RES,DUM) C IMPLICIT NONE C INTEGER NAF,NPF(*),NAP,NPP(*) INTEGER I,IDX,II,IX,IXP,IXX INTEGER J,JJ,K INTEGER NCY,NCZ,NPX,NPY,NPZ,NX,NY,NZ C REAL F(*),PSF(*),RES(*),DUM(*) C C init NX = NPF(1) IF (NAF.GE.2) THEN NY = NPF(2) ELSE NY = 1 ENDIF IF (NAF.GE.3) THEN NZ = NPF(3) ELSE NZ = 1 ENDIF NPX = NPP(1) IF (NAP.GE.2) THEN NPY = NPP(2) ELSE NPY = 1 ENDIF NCY = NPY / 2 + 1 IF (NAP.GE.3) THEN NPZ = NPP(3) ELSE NPZ = 1 ENDIF NCZ = NPZ / 2 + 1 C C for speed reasons use different code for 1, 2, 3-dimensions IF (NZ.EQ.1) GOTO 3000 C C 3-dim convolution IXX = 1 DO 2000, II=1,NZ C DO 1500, I=1,NY IX = IXX DO 400, K=1,NX RES(IX) = 0.0 IX = IX + 1 400 CONTINUE IXP = 1 C DO 1000, JJ=1,NPZ DO 800, J=1,NPY IDX = 1 + NX*(MAX(1,MIN(NY,J+I-NCY))-1 + + NY*(MAX(1,MIN(NZ,JJ+II-NCZ))-1) ) CALL CONVL1(F(IDX),NX,PSF(IXP),NPX,RES(IXX)) IXP = IXP + NPX 800 CONTINUE 1000 CONTINUE IXX = IXX + NX 1500 CONTINUE C 2000 CONTINUE C C we're done RETURN C 3000 IF (NY.EQ.1) GOTO 5000 C C handle 2-dim convolution IXX = 1 DO 4000, I=1,NY IX = IXX DO 3300, K = 1,NX RES(IX) = 0.0 IX = IX + 1 3300 CONTINUE IXP = 1 DO 3600 J=1,NPY IDX = 1 + NX*(MAX(1,MIN(NY,J+I-NCY))-1) CALL CONVL1(F(IDX),NX,PSF(IXP),NPX,RES(IXX)) IXP = IXP + NPX 3600 CONTINUE IXX = IXX + NX 4000 CONTINUE C C we're done RETURN C C handle 1-dim convolution 5000 IX = 1 DO 5500, K=1,NX RES(IX) = 0.0 IX = IX + 1 5500 CONTINUE CALL CONVL1(F,NX,PSF,NPX,RES) C C we're done RETURN C END SUBROUTINE CONVL1(F,NF,PSF,NP,C) C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION: C subroutine CONVL1 version 1.00 880916 C P. GROSBOL ESO - GARCHING C K. Banse (streamlined for MIDAS) C C.KEYWORDS: C convolution C C.PURPOSE: C convoles two one dimensional arrays with each other C C.ALGORITHM: C a one dimensional array is convoled with a one dimensional C point spread function. outside the limits of the function C array the first/last value is taken. C C.INPUT/OUTPUT: C the subroutine is called as CONVL1(F,NF,PSF,NP,C) C C F(NF) : function array (INPUT) C NF : dimension of function array (INPUT) C PSF(NP) : point spread function array (INPUT) C NP : dimension of point spread functon (INPUT) C C(NF) : convolved result array (OUTPUT) C C NOTE : the result is added to the output array 'C'. C C--------------------------------------------------------------- C IMPLICIT NONE C INTEGER NF,NP INTEGER I,K,NL1,NL2,NL3 INTEGER NOS,NR C REAL F(NF),C(NF),PSF(NP) REAL FACT,R C C initialize NOS = NP/2 NL1 = NOS + 1 NL2 = NF - NOS NL3 = NL2 + 1 C C if NP = 1, we only use middle section IF (NOS.EQ.0) THEN DO 1000, I=NL1,NL2 NR = I - NOS R = C(I) DO 400, K=1,NP R = R + PSF(K)*F(NR) NR = NR + 1 400 CONTINUE C(I) = R 1000 CONTINUE RETURN ENDIF C C for NP > 1, handle 1. section, where values before F(1) would be needed FACT = F(1) DO 2000, I=1,NOS R = C(I) DO 1400 K=1,1+NOS-I R = R + PSF(K)*FACT 1400 CONTINUE NR = 1 C DO 1600, K=2+NOS-I,NP R = R + PSF(K)*F(NR) NR = NR + 1 1600 CONTINUE C(I) = R 2000 CONTINUE C C in the middle no problems DO 3000 I=NL1,NL2 NR = I - NOS R = C(I) DO 2400, K=1,NP R = R + PSF(K)*F(NR) NR = NR + 1 2400 CONTINUE C(I) = R 3000 CONTINUE C C last section, values after F(NF) would be needed FACT = F(NF) DO 4000, I=NL3,NF NR = I - NOS R = C(I) DO 3300, K=1,NF+NOS-I R = R + PSF(K)*F(NR) NR = NR + 1 3300 CONTINUE C DO 3600, K=NF+NOS-I+1,NP R = R + PSF(K)*FACT 3600 CONTINUE C(I) = R 4000 CONTINUE C RETURN END REAL FUNCTION DECNVL(F,NAF,NPF,PSF,NAP,NPP,RES,BUF,DUM) C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION: C function DECNVL version 2.10 840116 C P. GROSBOL ESO - GARCHING C K. Banse (streamlined for MIDAS) C C.KEYWORDS: C DECONVOLUTION C C.PURPOSE: C deconvolve an image with its point spread function C C.ALGORITHM: C the iterative method of lucy is used for deconvolving upto C three dimensional images with their point spread function. C each call will performe one iterative step only. C REF.: LUCY,L.B. 1974 , A.J. VOL. 79, P. 745 C C.INPUT/OUTPUT: C call as DECNVL(F,NAF,NPF,PSF,NAP,NPP,RES,BUF,DUM) C C DECNVL : max. relative change applied C F : original image (INPUT) C NAF : dimension of F,RES and BUF (INPUT) C NPF() : no. of pixels in F,RES and BUF (INPUT) C PSF : point spread function (INPUT) C NAP : dimension of PSF (INPUT) C NPP : no. of pixels in PSF (INPUT) C RES : deconvolved image (INPUT/OUTPUT) C BUF : correction buffer (OUTPUT) C DUM() : dummy array (OUTPUT) C C--------------------------------------------------------------- C IMPLICIT NONE C INTEGER NAF,NAP,NPF(*),NPP(*) INTEGER I,IDX,II,IX,IXP,IXX INTEGER J,JJ,K,NCY,NCZ INTEGER NPX,NPY,NPZ,NX,NY,NZ C REAL F(*),RES(*),BUF(*),PSF(*),DUM(*) REAL ADUM,RELC C C init NX = NPF(1) IF (NAF.GE.2) THEN NY = NPF(2) ELSE NY = 1 ENDIF IF (NAF.GE.3) THEN NZ = NPF(3) ELSE NZ = 1 ENDIF NPX = NPP(1) IF (NAP.GE.2) THEN NPY = NPP(2) ELSE NPY = 1 ENDIF NCY = NPY / 2 + 1 IF (NAP.GE.3) THEN NPZ = NPP(3) ELSE NPZ = 1 ENDIF NCZ = NPZ / 2 + 1 RELC = 0.0 C C for speed reasons use different code for 1, 2, 3-dimensions IF (NZ.EQ.1) GOTO 3000 C C first convolution for 3-dim image IXX = 1 DO 1000, II=1,NZ DO 800, I=1,NY IX = IXX DO 200, K=1,NX BUF(IX) = 0.0 IX = IX + 1 200 CONTINUE IXP = 1 DO 600, JJ=1,NPZ DO 400, J=1,NPY IDX = 1 + NX*(MAX(1,MIN(NY,J+I-NCY))-1 + + NY*(MAX(1,MIN(NZ,JJ+II-NCZ))-1) ) CALL CONVL1(RES(IDX),NX,PSF(IXP),NPX,BUF(IXX)) IXP = IXP + NPX 400 CONTINUE 600 CONTINUE IX = IXX DO 700, K=1,NX BUF(IX) = F(IX) / BUF(IX) IX = IX + 1 700 CONTINUE IXX = IXX + NX 800 CONTINUE 1000 CONTINUE C C second convolution for 3-dim image DO 2400, II=1,NZ DO 2000, I=1,NY DO 1200, K=1,NX DUM(K) = 0.0 1200 CONTINUE IXP = 1 DO 1600, JJ=1,NPZ DO 1400, J=1,NPY IDX = 1 + NX*(MAX(1,MIN(NY,J+I-NCY))-1 + + NY*(MAX(1,MIN(NZ,JJ+II-NCZ))-1) ) CALL CONVL1(BUF(IDX),NX,PSF(IXP),NPX,DUM) IXP = IXP + NPX 1400 CONTINUE 1600 CONTINUE IX = IXX DO 1800, K=1,NX ADUM = ABS(1.-DUM(K)) IF (ADUM.GT.RELC) RELC = ADUM RES(IX) = DUM(K) * RES(IX) IX = IX + 1 1800 CONTINUE IXX = IXX + NX 2000 CONTINUE 2400 CONTINUE C C we're done, return DECNVL = RELC RETURN C C check y-dimension 3000 IF (NY.EQ.1) GOTO 6000 C C first convolution for 2-dim image IXX = 1 DO 4000, I=1,NY IX = IXX DO 3200, K=1,NX BUF(IX) = 0.0 IX = IX + 1 3200 CONTINUE IXP = 1 DO 3400, J=1,NPY IDX = 1 + NX*(MAX(1,MIN(NY,J+I-NCY))-1) CALL CONVL1(RES(IDX),NX,PSF(IXP),NPX,BUF(IXX)) IXP = IXP + NPX 3400 CONTINUE IX = IXX DO 3600, K=1,NX BUF(IX) = F(IX) / BUF(IX) IX = IX + 1 3600 CONTINUE IXX = IXX + NX 4000 CONTINUE C C second convolution for 2-dim image IXX = 1 DO 5000, I=1,NY DO 4200, K=1,NX DUM(K) = 0.0 4200 CONTINUE IXP = 1 DO 4400, J=1,NPY IDX = 1 + NX*(MAX(1,MIN(NY,J+I-NCY))-1) CALL CONVL1(BUF(IDX),NX,PSF(IXP),NPX,DUM) IXP = IXP + NPX 4400 CONTINUE IX = IXX DO 4600, K=1,NX ADUM = ABS(1.-DUM(K)) IF (ADUM.GT.RELC) RELC = ADUM RES(IX) = DUM(K) * RES(IX) IX = IX + 1 4600 CONTINUE IXX = IXX + NX 5000 CONTINUE C C we're done, return DECNVL = RELC RETURN C C first convolution for 1-dim image 6000 DO 6200, K=1,NX BUF(K) = 0.0 6200 CONTINUE CALL CONVL1(RES,NX,PSF,NPX,BUF) DO 6400, K=1,NX BUF(K) = F(K) / BUF(K) 6400 CONTINUE C C second convolution for 1-dim image DO 6600, K=1,NX DUM(K) = 0.0 6600 CONTINUE C CALL CONVL1(BUF,NX,PSF,NPX,DUM) DO 6800, K=1,NX ADUM = ABS(1.-DUM(K)) IF (ADUM.GT.RELC) RELC = ADUM RES(K) = DUM(K) * RES(K) 6800 CONTINUE C C we're done, return DECNVL = RELC RETURN C END