C @(#)backgrnd.for 17.1.1.1 (ES0-DMD) 01/25/02 17:19:19 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 C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C.IDENTIFICATION: BKGRND C.LANGUAGE: F77+ESOext C.AUTHOR: Capaccioli/Held C.KEYWORDS: polinomial fit C.PURPOSE: Calculate the coefficients of a 2D polinomial approximation C of the sky background. Optionally create frames with the C background subtracted image and the polinomial surface. C.ALGORITHM: Coefficients are calculated by an iterative least-squares C fitting of the pixels in the rebinned image frame. If a C predefined 'masking frame' is given, only not masked pixels C are taken in account. C.INPUT/OUTPUT: the following keywords are used: C IN_A/C/1/60 name of image data frame C IN_B/C/1/60 name of masking frame C OUT_A/C/1/60 name of output corrected image C OUT_B/C/1/60 name of output background surface C INPUTI/I/1/2 degree, maximum n. of iterations C INPUTR/R/1/1 clipping par. for the 1st iteration C INPUTR/R/2/1 clipping par. for following iter. C INPUTR/R/3/1 upper limit of acceptable skewness C INPUTC/C/1/8 name of the keyword storing coeff. C.VERSION 840517 C.VERSION: 871124 RHW ESO-FORTRAN Conversion C --------------------------------------------------------------------- PROGRAM BCKGRN IMPLICIT NONE C CHARACTER*60 FRAMEA,FRAMEB,FRAMEC,FRAMED CHARACTER*8 INPUTC CHARACTER CUNITA*64,IDENTA*72,CUNITB*64,IDENTB*72 CHARACTER HISTC*80,HISTD*80 CHARACTER MLINE1*40,MLINE2*40,LINE*80 INTEGER MADRID,KUN,KNUL,IMFA,IMFB,IMFD,IMFC INTEGER IAV,IPNTRA,IPNTRB,IPNTRC,IPNTRD,ISTAT INTEGER JER,MODE,NAXISA,NAXISB,NLOOP INTEGER INPUTI(2) INTEGER NPIXA(3),NPIXB(3),NPIXT(3) INTEGER ICONT,ICONTP C REAL AINCR,CONV REAL INPUTR(3) REAL ATEMP(128,128) REAL RKAPPA,REJ,RMEAN,SIGM1,SKEW1,RKURT DOUBLE PRECISION STEPA(3),STARTA(3),STEPB(3),STARTB(3) DOUBLE PRECISION STEPT(3),STARTT(3) DOUBLE PRECISION C(21),DSTORE(7) C LOGICAL FLMSK,FLGSB,FLGSU,CASEA,CASEB INCLUDE 'MID_INCLUDE:ST_DEF.INC/NOLIST' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC/NOLIST' C DATA MLINE1/' Iter Points Mean (10^6) Sigma'/ DATA MLINE2/' Skewness Kurtosis'/ DATA LINE/' '/ DATA CASEA,CASEB/.FALSE.,.FALSE./ DATA REJ/-99999./ ! minimum relative change DATA CONV/.002/ ! needed if no mask DATA NPIXB/1,1,1/ C C *** set up MIDAS environment + enable automatic error abort CALL STSPRO('BKGRND') CALL STKRDI('MODE',1,1,IAV,MODE,KUN,KNUL,ISTAT) C C *** get names of input frame + masking frame CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUN,KNUL,ISTAT) CALL STKRDC('IN_B',1,1,60,IAV,FRAMEB,KUN,KNUL,ISTAT) C C *** setup history update HISTC = 'background corrected frame obtained from '//FRAMEA HISTD = 'background pol. surface using dimensions from '// + FRAMEA C C *** get output frames CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,KUN,KNUL,ISTAT) CALL STKRDC('OUT_B',1,1,60,IAV,FRAMED,KUN,KNUL,ISTAT) C C *** read flag keyword INPUTC CALL STKRDC('INPUTC',1,1,12,IAV,INPUTC,KUN,KNUL,ISTAT) C C *** find out which case IF (INPUTC(1:1).EQ.'?') THEN ! read useful parameters CALL STKRDI('INPUTI',1,2,IAV,INPUTI,KUN,KNUL,ISTAT) CALL STKRDR('INPUTR',1,3,IAV,INPUTR,KUN,KNUL,ISTAT) CASEA = .TRUE. ELSE ! only create new images reading coeff. from a keyword CASEB = .TRUE. END IF C C *** get input frame and map it CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, 2 2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA, 3 IPNTRA,IMFA,ISTAT) IF (NAXISA.GT.2) THEN CALL STTPUT('*** FATAL: maximum 2 dim. currently supported...', + ISTAT) CALL STSEPI END IF C C *** get masking frame (if exists) and map it IF (FRAMEB(1:1).NE.'+') THEN FLMSK = .TRUE. CALL STIGET(FRAMEB,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE, 2 2,NAXISB,NPIXB,STARTB,STEPB,IDENTB,CUNITB, 3 IPNTRB,IMFB,ISTAT) ELSE FLMSK = .FALSE. END IF C C *** get output frame for corrected image, map it and set descriptors IF (FRAMEC(1:1).NE.'+') THEN FLGSB = .TRUE. CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, 2 NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA, 3 IPNTRC,IMFC,ISTAT) CALL STDCOP(IMFA,IMFC,4,'HISTORY',ISTAT) CALL STDCOP(IMFA,IMFC,4,'LHCUTS',ISTAT) C CALL STDWRC(IMFC,'HISTORY',1,HISTC,-1,80,ISTAT) ELSE FLGSB = .FALSE. END IF C C *** get output frame for polinomial surface , map it and set descriptors IF (FRAMED(1:1).NE.'+') THEN FLGSU = .TRUE. CALL STIPUT(FRAMED,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, 2 NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA, 3 IPNTRD,IMFD,ISTAT) CALL STDCOP(IMFA,IMFD,4,'HISTORY',ISTAT) CALL STDCOP(IMFA,IMFD,4,'LHCUTS',ISTAT) C CALL STDWRC(IMFD,'HISTORY',1,HISTD,-1,80,ISTAT) ELSE FLGSU = .FALSE. END IF C C *** check parameters values IF (CASEA) THEN IF (INPUTI(1).LE.0 .OR. INPUTI(1).GT.5) THEN CALL STTPUT('*** FATAL: Invalid degree: should be 1 to 5', + ISTAT) CALL STSEPI END IF IF (INPUTI(2).LE.0) THEN CALL STTPUT('*** FATAL: Invalid number of iterations: '// + 'should be greater than zero',ISTAT) CALL STSEPI END IF IF (INPUTI(2).GT.20) THEN INPUTI(2) = 20 CALL STTPUT('*** WARNING: Maximum number of iterations '// + 'set to 20 ',ISTAT) END IF IF (INPUTR(1).LE.0 .OR. INPUTR(2).LE.0) THEN CALL STTPUT('*** FATAL: Invalid kappa: should be '// + 'greater than zero',ISTAT) CALL STSEPI END IF IF (INPUTR(3).LE.0) THEN CALL STTPUT('*** FATAL: Invalid skewness: should be '// + 'greater than zero',ISTAT) CALL STSEPI END IF END IF C C *** if case b jump to section creating output frames IF (CASEB) GO TO 20 C C *** rebin input image to run faster CALL SQZTMP(MADRID(IPNTRA),ATEMP,NPIXA(1),NPIXA(2),STARTA(1), + STARTA(2),STEPA(1),STEPA(2),NPIXT(1),NPIXT(2), + STARTT(1),STARTT(2),STEPT(1),STEPT(2)) C C *** write scaling parameters and degree in keyword INPUTC DSTORE(1) = NPIXT(1) DSTORE(2) = NPIXT(2) DSTORE(3) = STARTT(1) DSTORE(4) = STARTT(2) DSTORE(5) = STEPT(1) DSTORE(6) = STEPT(2) DSTORE(7) = INPUTI(1) CALL STKWRD('BG_COEF',DSTORE,22,7,KUN,ISTAT) C C *** main loop : iterations RKAPPA = INPUTR(1) DO 10 NLOOP = 1,INPUTI(2) IF (NLOOP.NE.1) RKAPPA = INPUTR(2) C C *** calculate coefficients for this loop and param. of residuals CALL FEET(MADRID(IPNTRB),ATEMP,INPUTI(1),FLMSK,REJ,NPIXB(1), + NPIXB(2),STARTB(1),STARTB(2),STEPB(1),STEPB(2), + NPIXT(1),NPIXT(2),STARTT(1),STARTT(2),STEPT(1), + STEPT(2),ICONT,RMEAN,SIGM1,SKEW1,RKURT,C,JER) C C *** check interpolation failure IF (JER.NE.0) THEN CALL STTPUT('*** FATAL: Interpolation failed',ISTAT) CALL STSEPI END IF C C *** write coefficients in keyword INPUTC CALL STKWRD('BG_COEF',C,1,21,KUN,ISTAT) C C *** write significant parameters on terminal IF (NLOOP.EQ.1) THEN CALL STTPUT(MLINE1//MLINE2,ISTAT) ENDIF WRITE(LINE(1:16),'(2I8)') NLOOP,ICONT CALL VFL(RMEAN*1.E6,LINE(19:26)) CALL VFL(SIGM1,LINE(29:36)) CALL VFL(SKEW1,LINE(39:46)) CALL VFL(RKURT,LINE(49:56)) CALL STTPUT(LINE,ISTAT) C C *** if skewness low enough , stop iteration IF (ABS(SKEW1).LE.INPUTR(3)) THEN CALL STSEPI ENDIF C C *** also if no significant change occurred IF (NLOOP.NE.1) THEN AINCR = FLOAT(ICONTP-ICONT)/ICONTP IF (AINCR.LT.CONV) THEN CALL STTPUT('*** WARNING: iteration stopped: not '// + ' enough pixels were rejected',ISTAT) GO TO 20 END IF END IF ICONTP = ICONT ! save n. of processed points C C *** if last iteration , skip mask update IF (NLOOP.EQ.INPUTI(2)) GO TO 20 C C *** if not enough pixels , stop iteration IF (ICONT.LT.50) THEN CALL STTPUT('*** WARNING: Iteration stopped: '// 2 'not enough points ',ISTAT) GO TO 20 END IF C C *** now update masking frame and temporary buffer for next iteration CALL UPDAT(MADRID(IPNTRB),ATEMP,FLMSK,REJ,NPIXB(1),NPIXB(2), + STARTB(1),STARTB(2),STEPB(1),STEPB(2),NPIXT(1), + NPIXT(2),STARTT(1),STARTT(2),STEPT(1),STEPT(2), + INPUTI(1),RKAPPA,SIGM1,C) C 10 CONTINUE ! end main loop C 20 CONTINUE ! this is the skip out point C C *** read coefficients from keyword (case b) IF (CASEB) THEN CALL STKRDD('BG_COEF',1,21,C,IAV,KUN,ISTAT) CALL STKRDD('BG_COEF',22,7,DSTORE,IAV,KUN,ISTAT) NPIXT(1) = DSTORE(1) NPIXT(2) = DSTORE(2) STARTT(1) = DSTORE(3) STARTT(2) = DSTORE(4) STEPT(1) = DSTORE(5) STEPT(2) = DSTORE(6) INPUTI(1) = DSTORE(7) END IF C C *** create back.corrected and/or pol. surf. image frame (testing flags) IF (.NOT.(FLGSB.OR.FLGSU)) THEN CALL STSEPI ELSE CALL SBTRCT(MADRID(IPNTRA),MADRID(IPNTRC),MADRID(IPNTRD), + NPIXA(1),NPIXA(2),STARTA(1),STARTA(2),STEPA(1), + STEPA(2),NPIXT(1),NPIXT(2),STARTT(1),STARTT(2), + STEPT(1),STEPT(2),C,INPUTI(1),FLGSB,FLGSU) CALL STSEPI ENDIF END SUBROUTINE SBTRCT(FRAMEA,FRAMEC,FRAMED,NIX,NIY,DXSTR,DYSTR, + DXSTP,DYSTP,NSX,NSY,SXSTR,SYSTR,SXSTP,SYSTP,C, + IDEG,FLGSB,FLGSU) C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++: C.PURPOSE: Create new image frames with the background subtracted image and/or C the polinomial surface ,testing option flags FLGSB and FLGSU: C.INPUT: FRAMEA : INPUT_A is the input image C FLGSB : flag = .TRUE. if you want backgrond subtraction C FLGSU : flag = .TRUE. if you want a background surface image C.OUTPUT: FRAMEC : OUTPUT_C is the background subtracted image C FRAMED : OUTPUT_D is the polinomial surface C ------------------------------------------------------------------ IMPLICIT NONE REAL FRAMEA(1) REAL FRAMEC(1) REAL FRAMED(1) INTEGER NIX INTEGER NIY DOUBLE PRECISION DXSTR DOUBLE PRECISION DYSTR DOUBLE PRECISION DXSTP DOUBLE PRECISION DYSTP INTEGER NSX INTEGER NSY DOUBLE PRECISION SXSTR DOUBLE PRECISION SYSTR DOUBLE PRECISION SXSTP DOUBLE PRECISION SYSTP DOUBLE PRECISION C(21) INTEGER IDEG LOGICAL FLGSB LOGICAL FLGSU C INTEGER I, IX, IY, IF INTEGER J REAL DCX, DCY REAL POLIN REAL Y1, YSRV REAL X1, XSRV DOUBLE PRECISION CC(6) DOUBLE PRECISION XD,YD C C *** initialize variables C DCX = (NSX-1.)/2. DCY = (NSY-1.)/2. C C *** fill image(S) C IF = 1 DO 20 J = 1,NIY IY = J - 1 Y1 = DYSTR + DYSTP*FLOAT(IY) C C *** find pixel coord. in temporary array used by fitting routine C YSRV = (Y1-SYSTR)/SYSTP YD = DCY - YSRV C C *** initialize calculation of coefficients C CALL COEFY(YD,C,CC,IDEG) C DO 10 I = IF,IF + NIX - 1 IX = I - IF X1 = DXSTR + DXSTP*FLOAT(IX) XSRV = (X1-SXSTR)/SXSTP XD = XSRV - DCX IF (FLGSB) THEN FRAMEC(I) = FRAMEA(I) - POLIN(XD,CC,IDEG) IF (FLGSU) THEN FRAMED(I) = POLIN(XD,CC,IDEG) ENDIF ELSE IF (FLGSU) THEN FRAMED(I) = POLIN(XD,CC,IDEG) ENDIF END IF 10 CONTINUE C IF = IF + NIX 20 CONTINUE C RETURN END