* Last processed by NICE on 12-Jun-2000 15:51:00
* Customized for :  IEEE, LINUX, UNIX, MOTIF, F77
      PROGRAM BACKGROUND_MAIN
C----------------------------------------------------------------------
C Compute a large scale background
C----------------------------------------------------------------------
      INCLUDE 'inc:xpar.inc'
      INCLUDE 'inc:ypar.inc'
      INCLUDE 'inc:format.inc'
      INCLUDE 'inc:errcod.inc'
      INCLUDE 'inc:memory.inc'
*
      CHARACTER*80 NAMEX,NAMEY
      CHARACTER*4 METHOD
      INTEGER MODE
      LOGICAL ERROR,SLOW,FOURIER
      INTEGER N, I, SCALE, NW, NT, LENC, SIC_GETVM
      INTEGER NHISTO,NBACK,IER
      REAL PMIN,PMAX,PINC,X1,X2
      INTEGER*4 IPX, IPY, POINTER, ADDR, X,Y,Z,W,IW,HISTO
      INTEGER*4 BLC(4), TRC(4)
      INTEGER NX,NY
      INTEGER MX,MY
      REAL CVI(6),CVO(6)
      COMMON/CONV/CVI,CVO
*
      DATA BLC /4*0/, TRC /4*0/
*
      SLOW=.FALSE.
      CALL GILDAS_OPEN
      CALL GILDAS_CHAR('Y_NAME$',NAMEY)
      CALL GILDAS_CHAR('X_NAME$',NAMEX)
      CALL GILDAS_REAL('X_MIN$',X1,1)
      CALL GILDAS_REAL('X_MAX$',X2,1)
      CALL GILDAS_INTE('NBINS$',NHISTO,1)
      CALL GILDAS_INTE('WIDTH$',NW,1)
      CALL GILDAS_INTE('SCALE$',SCALE,1)
      CALL GILDAS_INTE('N_MIN$',NT,1)
      CALL GILDAS_CHAR('METHOD$',METHOD)
      CALL GILDAS_CLOSE
*
      N = LENC(NAMEY)
      CALL SIC_PARSEF(NAMEY(1:N),Y_FILE,' ','.gdf')
      CALL GDF_GEIS (Y_ISLO,ERROR)
      IF (.NOT.ERROR) CALL
     $GDF_REIS (Y_ISLO,Y_TYPE,Y_FILE,Y_FORM,Y_SIZE,ERROR)
      IF (ERROR) THEN
         WRITE(6,*) 'F-BACKGROUND,  Cannot read input file'
         GOTO 100
      ENDIF
      IF (Y_FORM.NE.FMT_R4) THEN
         WRITE(6,*) 'E-BACKGROUND,  Work only on Real images'
         GOTO 100
      ENDIF
      CALL GDF_READY (Y_ISLO,ERROR)
*
      CALL GDF_CHYX
*
      NX=X_DIM(1)
      NY=X_DIM(2)
      PMIN = MIN(X1,X2)
      PMAX = MAX(X1,X2)
      PINC = (PMAX-PMIN)/NHISTO
*
* Consistency checks.
      NT = MAX(NT,3*NHISTO)
      NT = MIN(NT,3*NW**2/4)             ! PI*D^2/4
      MODE=1
      SLOW=.FALSE.
      FOURIER=.FALSE.
      IF (METHOD.EQ.'ZERO') THEN
         MODE=0
      ELSEIF (METHOD.EQ.'FAST') THEN
         MODE=1
      ELSEIF (METHOD.EQ.'LINE') THEN
         MODE=1
      ELSEIF (METHOD.EQ.'PARA') THEN
         MODE=2
      ELSEIF (METHOD.EQ.'GAUS') THEN
         MODE=3
      ELSEIF (METHOD.EQ.'FOUR') THEN
         FOURIER = .TRUE.
         MODE=-1
      ELSEIF (METHOD.EQ.'SLOW') THEN
         SLOW=.TRUE.
         MODE=0
      ENDIF
*
* Get some work space
      IER = SIC_GETVM(NHISTO,ADDR)
      IF (IER.NE.1) GOTO 99
      HISTO = POINTER(ADDR,MEMORY)
      MX=INT(FLOAT(NX-NW)/FLOAT(SCALE))+1
      MY=INT(FLOAT(NY-NW)/FLOAT(SCALE))+1
      DO I=1,6
         CVI(I)=1.0
      ENDDO
      CVO(1)=1.0
      CVO(2)=FLOAT(NW+1)*0.5
      CVO(3)=FLOAT(SCALE)
      CVO(4)=1.0
      CVO(5)=FLOAT(NW+1)*0.5
      CVO(6)=FLOAT(SCALE)
      NBACK = MX*MY
      IF (SLOW) THEN
         WRITE(6,*) 'Using slow method...'
         IER = SIC_GETVM(3*NBACK,ADDR)
         IF (IER.NE.1) GOTO 99
         X = POINTER(ADDR,MEMORY)
         Y = X + NBACK
         Z = X + 2*NBACK
         IER = SIC_GETVM(35*NBACK,ADDR)
         IF (IER.NE.1) GOTO 99
         IW = POINTER(ADDR,MEMORY)
         IER = SIC_GETVM(13*NBACK,ADDR)
         IF (IER.NE.1) GOTO 99
         W = POINTER(ADDR,MEMORY)
      ELSEIF (FOURIER) THEN
         X_DIM(1)=MX
         X_DIM(2)=MY
         X_SIZE = X_DIM(1)*X_DIM(2)*X_DIM(3)*X_DIM(4)
      ELSE
         IER = SIC_GETVM(3*NBACK,ADDR)   ! z(i,j),xgrav(i,j),ygrav(i,j)
         IF (IER.NE.1) GOTO 99
         X = POINTER(ADDR,MEMORY)
      ENDIF
      N = LENC(NAMEX)
      CALL SIC_PARSEF(NAMEX(1:N),X_FILE,' ','.gdf')
      CALL GDF_GEIS (X_ISLO,ERROR)
      IF (ERROR) THEN
         WRITE(6,*) 'F-BACKGROUND,  Cannot create output image'
         GOTO 100
      ENDIF
      CALL GDF_WRITX (X_ISLO,ERROR)
*
* Create output image
      CALL GDF_CRIS (X_ISLO,X_TYPE,X_FILE,X_FORM,X_SIZE,ERROR)
      IF (ERROR) THEN
         WRITE(6,*) 'F-BACKGROUND,  Cannot create output image'
         GOTO 100
      ENDIF
*
* Load the data into memory
      CALL GDF_GEMS (X_MSLO,X_ISLO,BLC,TRC,X_ADDR,X_FORM,ERROR)
      IF (ERROR) GOTO 100
      IPX = POINTER(X_ADDR,MEMORY)
      CALL GDF_GEMS (Y_MSLO,Y_ISLO,BLC,TRC,Y_ADDR,Y_FORM,ERROR)
      IF (ERROR) GOTO 100
      IPY = POINTER(Y_ADDR,MEMORY)
*
* Do the real work
      IF (SLOW) THEN
         CALL ZODI011 (MEMORY(IPY),      !Input image
     $   NX,NY,                          !Size
     $   X_DIM(3),
     $   NT,                             !Minimum number of data points
     $   SCALE,                          !Scale
     $   NW,                             !Width
     $   X_BVAL,                         !Blanking value
     $   X_EVAL,
     $   PMIN,PMAX,PINC,
     $   MEMORY(HISTO),
     $   MEMORY(X),
     $   MEMORY(Y),
     $   MEMORY(Z),
     $   MEMORY(W),
     $   MEMORY(IW),
     $   MEMORY(IPX))
      ELSE
* Define the background map
         IF (X_BLAN.EQ.0) THEN
            X_BVAL=PMIN-1.0
            X_EVAL=0.9
         ENDIF
         IF (FOURIER) X = IPX
         CALL ZODI001 (MEMORY(IPY),      !Input image
     $   NX,NY,                          !Size
     $   Y_DIM(3),
     $   NT,                             !Minimum number of data points
     $   SCALE,                          !Scale
     $   NW,                             !Width
     $   X_BVAL,                         !Blanking value
     $   X_EVAL,
     $   PMIN,PMAX,PINC,
     $   MEMORY(HISTO),
     $   MEMORY(X),
     $   MEMORY(IPX),
     $   MODE)
      ENDIF
* Update disk files and free memory
      CALL GDF_FRMS(X_MSLO,ERROR)
      CALL GDF_FRMS(Y_MSLO,ERROR)
*
      STOP 'S-BACKGROUND,  Successful completion'
*
99    WRITE(6,*) 'F-BACKGROUND,  Memory allocation failure'
100   CALL SYSEXI (FATALE)
      END
*<FF>
      SUBROUTINE ZODI001(A,NX,NY,NZ,NT,SCALE,NW,BVAL,EVAL,
     $PMIN,PMAX,PINC,HISTO,INTER,B,MODE)
C----------------------------------------------------------------------
C Compute a large scale background, FAST VERSION
C
C Arguments:
C	A	R*4	Input data cube				Input
C	NX	I	First dimension of A			Input
C	NY	I	Second dimension of A			Input
C	NZ	I	Third dimension of A			Input
C	NT	I	Minimum valid number of points 		Input
C	SCALE	I	Scale factor for intermediate grid	Input
C	NW	I	Diameter of smoothing circle		Input
C	BVAL	R*4	Blanking value				Input
C	EVAL	R*4	Tolerance on blanking value		Input
C	PMIN	R*4	Minimum accepted value			Input
C	PMAX	R*4	Maximum accepted value			Input
C	PINC	R*4	Histogram step				Input
C	HISTO	R*4	Histogram array				Scratch
C	INTER	R*4	Intermediate image			Scratch
C	B	R*4	Output background image			Output
C       MODE    I       0=nointerpol, 1 or 2 interpol same orderInput
C----------------------------------------------------------------------
*
      INTEGER NX,NY,NZ,MX,MY,SCALE,NW,IER,NT,MODE
      REAL*4 A(NX,NY,NZ),B(NX,NY,NZ),INTER(*)
      REAL*4 HISTO(*),BVAL,PMIN,PMAX,PINC,EVAL
*
      INTEGER NH,I
*
      MX=INT(FLOAT(NX-NW)/FLOAT(SCALE))+1
      MY=INT(FLOAT(NY-NW)/FLOAT(SCALE))+1
* Define the Histogram
      NH = NINT((PMAX-PMIN)/PINC)
*
* Compute the background image
      DO I=1,NZ                          ! Loop on planes
         CALL ZODI002(A(1,1,I),NX,NY,SCALE,NW,PMIN,
     $   PMAX,PINC,HISTO,NH,NT,INTER,MX,MY,BVAL)
*
* Resample the image to the full Grid
         IF (MODE.NE.-1) CALL ZODI003(INTER,MX,MY,
     $   B(1,1,I),NX,NY,BVAL,EVAL,
     $   SCALE,NW,MODE)
*
* Subtract the background from the image
*         CALL ZODI004 (A,B,NX,NY,BVAL,EVAL)
      ENDDO
      END
*<FF>
      SUBROUTINE ZODI002(A,NX,NY,NS,NW,PMIN,PMAX,PINC,H,NH,NT,
     $INTER,MX,MY,BVAL)
C----------------------------------------------------------------------
C Estimate large scale background on reduced grid, FAST VERSION
C Arguments:
C	A	R*4	Input data cube				Input
C	NX	I	First dimension of A			Input
C	NY	I	Second dimension of A			Input
C	NS	I	Scale factor for intermediate grid	Input
C	NW	I	Diameter of smoothing circle		Input
C	PMIN	R*4	Minimum accepted value			Input
C	PMAX	R*4	Maximum accepted value			Input
C	PINC	R*4	Histogram step				Input
C	H	R*4	Histogram array				Scratch
C	NH	I	Size of histogram			Input
C	NT	I	Minimum valid number of points 		Input
C     INTER	R*4	Intermediate image		        Output
C	MX	I	First dimension of INTER		Input
C	MY	I	Second dimension of INTER		Input
C	BVAL	R*4	Blanking value				Input
C----------------------------------------------------------------------
      INTEGER NX,NY,NH,M,NW,NS,NT,MX,MY
      REAL*4 A(NX,NY),H(NH),PMIN,PMAX,PINC,INTER(MX,MY,3),BVAL
      REAL*4 JA,IA
*
      REAL R,GRAV(2)
      INTEGER I,J,NF,K,IER
*
      REAL CVI(6),CVO(6)
      COMMON/CONV/CVI,CVO
*
      K = 0
      DO J = 1,MY
         JA = (J-CVO(4))*CVO(6)+CVO(5)
         DO I = 1,MX
            IA = (I-CVO(1))*CVO(3)+CVO(2)
            CALL ZODI102(A,NX,NY,IA,JA,PMIN,PMAX,PINC,
     $      H,NH,NW,NF)
            IF (NF.GE.NT) THEN
               CALL ZODI122(H,NH,R,PMIN,PINC,IER)
               IF (IER.EQ.0) THEN
                  CALL ZODI123(A,NX,NY,IA,JA,NW,PINC,R,GRAV,IER)
                  IF (IER.EQ.0) THEN
                     K    = K+1
                     INTER(I,J,1) = R
                     INTER(I,J,2) = GRAV(1)
                     INTER(I,J,3) = GRAV(2)
                  ENDIF
               ELSE
                  INTER(I,J,1) = BVAL
               ENDIF
            ENDIF
         ENDDO
      ENDDO
      M = K
      WRITE(6,*) M,' valid background values found'
      END
*<FF>
      SUBROUTINE ZODI003(A,MX,MY,B,NX,NY,BVAL,EVAL,NS,NW,MODE)
C----------------------------------------------------------------------
C	Regrid the intermediate background map on the same scale as
C       the original map, zero,1 or 2-order linear interpolation
C Arguments
C       A	R*4(*)	Input map                               Input
C	MX	I	First dimension of A			Input
C	MY	I	Second dimension of A			Input
C	B	R*4(*)	Output map                              Output
C	NX	I	First dimension of B			Input
C	NY	I	Second dimension of B			Input
C	BVAL	R*4	Blanking value                          Input
C	EVAL	R*4	Tolerance on blanking value		Input
C	NS	I	Scale factor for intermediate grid	Input
C	NW	I	Diameter of smoothing circle		Input
C       MODE    I       0=nointerpol, 1 or 2 interpol same orderInput
C----------------------------------------------------------------------
      INTEGER NX,NY,MX,MY,NS,MODE
      REAL*4 A(MX,MY,3), B(NX,NY), BVAL, EVAL
*
      INTEGER IB,JB,I,J,IMIN,IMAX,JMIN,JMAX,NW
      REAL XA,YA,XR,YR,APLUS,AZERO,AMOIN
      INTEGER IA,JA
*
      REAL X,Y,Z,X0,Y0,DX,DY,TMP,VALUE,WEIGHT,DELTA,DIST
      INTEGER IBB,JBB,ISUP
      REAL MYEXP(100)
*
      REAL CVI(6),CVO(6)
      COMMON/CONV/CVI,CVO
*
* Interpolate,Take Blanking into account, version #1
* Loop over Output data points
* blank output map
      DO JB=1,NY
         DO IB=1,NX
            B(IB,JB) = BVAL
         ENDDO
      ENDDO
*
* Gauss mode
      IF (MODE.EQ.3) THEN
         DELTA=FLOAT(NW)*FLOAT(NW)/(4*ALOG(2.0))
         DO I=1,100
            X=(I-1)*0.03*FLOAT(NW)
            X=X*X
            MYEXP(I)=EXP(-1.0*X/DELTA)
            WRITE (6,*) (I-1)*0.03,MYEXP(I)
         ENDDO
         ISUP=INT(FLOAT(NW)/CVO(3))
         DO JB=1,NY
            Y0=FLOAT(JB)
            TMP  = (JB-CVO(5))/CVO(6)+CVO(4)
            J  = NINT(TMP)
            JMIN = MAX(1,J-ISUP)         !support de la fonction Gauss dans A
            JMAX = MIN(MY,J+ISUP)
            DO IB=1,NX
               X0=FLOAT(IB)
               TMP  = (IB-CVO(2))/CVO(3)+CVO(1)
               I  = NINT(TMP)
               IMIN = MAX(1,I-ISUP)
               IMAX = MIN(MX,I+ISUP)
               VALUE=0.0
               WEIGHT=0.0
               DO JBB=JMIN,JMAX
                  DO IBB=IMIN,IMAX
                     Z=A(IBB,JBB,1)
                     IF (Z.NE.BVAL) THEN
                        X=A(IBB,JBB,2)
                        Y=A(IBB,JBB,3)
                        DX=(X-X0)/FLOAT(NW)
                        DY=(Y-Y0)/FLOAT(NW)
                        DIST=DX*DX+DY*DY
                        TMP=MYEXP(INT(100.0*DIST)+1)
                        VALUE=VALUE+Z*TMP
                        WEIGHT=WEIGHT+TMP
                     ENDIF
                  ENDDO
               ENDDO
               IF (WEIGHT.GT.0.0) THEN
                  B(IB,JB)=VALUE/WEIGHT
               ELSE
                  B(IB,JB)=BVAL
               ENDIF
            ENDDO
         ENDDO
      ELSEIF (MODE.EQ.0) THEN
         TMP  = (1-CVO(1))*CVO(3)+CVO(2)
         IMIN = NINT(TMP)
         TMP  = (MX-CVO(1))*CVO(3)+CVO(2)
         IMAX = NINT(TMP)
         TMP  = (1-CVO(4))*CVO(6)+CVO(5)
         JMIN = NINT(TMP)
         TMP  = (MY-CVO(4))*CVO(6)+CVO(5)
         JMAX = NINT(TMP)
         DO JB = JMIN,JMAX
            DO IB = IMIN,IMAX
               TMP  = (IB-CVO(2))/CVO(3)+CVO(1)
               I  = NINT(TMP)
               TMP  = (JB-CVO(5))/CVO(6)+CVO(4)
               J  = NINT(TMP)
               B(IB,JB) = A(I,J,1)
            ENDDO
         ENDDO
      ELSEIF (MODE.EQ.1) THEN
         TMP  = (1-CVO(1))*CVO(3)+CVO(2)
         IMIN = NINT(TMP)
         TMP  = (MX-1-CVO(1))*CVO(3)+CVO(2)
         IMAX = NINT(TMP)
         TMP  = (1-CVO(4))*CVO(6)+CVO(5)
         JMIN = NINT(TMP)
         TMP  = (MY-1-CVO(4))*CVO(6)+CVO(5)
         JMAX = NINT(TMP)
         DO JB = JMIN,JMAX
            DO IB = IMIN,IMAX
               TMP  = (IB-CVO(2))/CVO(3)+CVO(1)
               I  = NINT(TMP)
               TMP  = (JB-CVO(5))/CVO(6)+CVO(4)
               J  = NINT(TMP)
               XA = (I-CVO(1))*CVO(3)+CVO(2)
               YA = (J-CVO(4))*CVO(6)+CVO(5)
               IA = NINT(XA)
               JA = NINT(YA)
               XR = (FLOAT(IB)-FLOAT(IA))/CVO(3)+0.5
               YR = (FLOAT(JB)-FLOAT(JA))/CVO(6)+0.5
               IF (A(I,J,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I+1,J,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I,J+1,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I+1,J+1,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSE
                  B(IB,JB) = (1-XR)*(1-YR)*A(I,J,1)+
     $            XR*(1-YR)*A(I+1,J,1)+
     $            XR*YR*A(I+1,J+1,1)+
     $            (1-XR)*YR*A(I,J+1,1)
               ENDIF
            ENDDO
         ENDDO
      ELSEIF (MODE.EQ.2) THEN
         TMP  = (2-CVO(1))*CVO(3)+CVO(2)
         IMIN = NINT(TMP)
         TMP  = (MX-1-CVO(1))*CVO(3)+CVO(2)
         IMAX = NINT(TMP)
         TMP  = (2-CVO(4))*CVO(6)+CVO(5)
         JMIN = NINT(TMP)
         TMP  = (MY-1-CVO(4))*CVO(6)+CVO(5)
         JMAX = NINT(TMP)
         DO JB = JMIN,JMAX
            DO IB = IMIN,IMAX
               TMP  = (IB-CVO(2))/CVO(3)+CVO(1)
               I  = NINT(TMP)
               TMP  = (JB-CVO(5))/CVO(6)+CVO(4)
               J  = NINT(TMP)
               XA = (I-CVO(1))*CVO(3)+CVO(2)
               YA = (J-CVO(4))*CVO(6)+CVO(5)
               IA = NINT(XA)
               JA = NINT(YA)
               XR = (FLOAT(IB)-FLOAT(IA))/CVO(3)+0.5
               YR = (FLOAT(JB)-FLOAT(JA))/CVO(6)+0.5
               IF (A(I-1,J-1,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I,J-1,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I+1,J-1,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I-1,J,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I,J,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I+1,J,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I-1,J+1,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I,J+1,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSEIF (A(I+1,J+1,1).EQ.BVAL) THEN
                  B(IB,JB) = BVAL
               ELSE
*
* Interpolate (XR or YR first, does not matter in this case)
                  APLUS = XR*(XR-1)*A(I-1,J+1,1)*0.5+
     $            XR*(XR+1)*A(I+1,J+1,1)*0.5-
     $            (XR+1)*(XR-1)*A(I,J+1,1)
                  AZERO = XR*(XR-1)*A(I-1,J,1)*0.5+
     $            XR*(XR+1)*A(I+1,J,1)*0.5-
     $            (XR+1)*(XR-1)*A(I,J,1)
                  AMOIN = XR*(XR-1)*A(I-1,J-1,1)*0.5+
     $            XR*(XR+1)*A(I+1,J-1,1)*0.5-
     $            (XR+1)*(XR-1)*A(I,J-1,1)
*     Then YR (or XR)
                  B(IB,JB) = YR*(YR-1)*AMOIN*0.5+
     $            YR*(YR+1)*APLUS*0.5-
     $            (YR+1)*(YR-1)*AZERO
               ENDIF
            ENDDO
         ENDDO
      ENDIF
      END
*<FF>
      SUBROUTINE ZODI011(A,NX,NY,NZ,NT,SCALE,NW,BVAL,EVAL,
     $PMIN,PMAX,PINC,HISTO,X,Y,Z,W,IW,B)
C----------------------------------------------------------------------
C Compute a large scale background, SLOW METHOD
C
C Arguments:
C	A	R*4	Input data cube				Input
C	NX	I	First dimension of A			Input
C	NY	I	Second dimension of A			Input
C	NZ	I	Third dimension of A			Input
C	NT	I	Minimum valid number of points 		Input
C	SCALE	I	Scale factor for intermediate grid	Input
C	NW	I	Diameter of smoothing circle		Input
C	BVAL	R*4	Blanking value				Input
C	EVAL	R*4	Tolerance on blanking value		Input
C	PMIN	R*4	Minimum accepted value			Input
C	PMAX	R*4	Maximum accepted value			Input
C	PINC	R*4	Histogram step				Input
C	HISTO	R*4	Histogram array				Scratch
C	X	R*4	Work space				Scratch
C	Y	R*4	Work space				Scratch
C	Z	R*4	Work space				Scratch
C	W	R*4	Work space				Scratch
C	IW	R*4	Work space				Scratch
C	B	R*4	Output background image			Output
C----------------------------------------------------------------------
*
      INTEGER NX,NY,NZ,SCALE,NW,IER,NT
      REAL*4 A(NX,NY,NZ),B(NX,NY,NZ),X(*),Y(*),Z(*),W(*),IW(*)
      REAL*4 HISTO(*),BVAL,PMIN,PMAX,PINC,EVAL
*
      INTEGER MX,MY,M,NH,I
*
* Define the background map
      MX=INT(FLOAT(NX-NW)/FLOAT(SCALE))+1
      MY=INT(FLOAT(NY-NW)/FLOAT(SCALE))+1
      NH = NINT((PMAX-PMIN)/PINC)
*
* Compute the background image
      M = MX*MY
      DO I=1,NZ                          ! Loop on planes
         CALL ZODI012(A(1,1,I),NX,NY,SCALE,NW,PMIN,
     $   PMAX,PINC,HISTO,NH,NT,
     $   X,Y,Z,MX,MY,M)
*
* Resample the image to the full Grid
         CALL ZODI013(M,X,Y,Z,W,IW,B(1,1,I),NX,NY,BVAL,IER)
         IF (IER.NE.0) RETURN
*
* Subtract the background from the image
*         CALL ZODI004 (A,B,NX,NY,BVAL,EVAL)
      ENDDO
      END
*<FF>
      SUBROUTINE ZODI012(A,NX,NY,NS,NW,PMIN,PMAX,PINC,H,NH,NT,
     $X,Y,Z,MX,MY,M)
C----------------------------------------------------------------------
C Estimate large scale background on reduced grid
C
C Arguments:
C	A	R*4	Input data cube				Input
C	NX	I	First dimension of A			Input
C	NY	I	Second dimension of A			Input
C	NS	I	Scale factor for intermediate grid	Input
C	NW	I	Diameter of smoothing circle		Input
C	PMIN	R*4	Minimum accepted value			Input
C	PMAX	R*4	Maximum accepted value			Input
C	PINC	R*4	Histogram step				Input
C	H	R*4	Histogram array				Scratch
C	NH	I	Size of histogram			Input
C	NT	I	Minimum valid number of points 		Input
C	BVAL	R*4	Blanking value				Input
C	EVAL	R*4	Tolerance on blanking value		Input
C	X	R*4	X coordinate for reduced grid		Output
C	Y	R*4	Y coordinate for reduced grid		Output
C	Z	R*4	Background value on reduced grid	Output
C	MX	I	First dimension of submap		Input
C	MY	I	Second dimension of submap		Input
C	M	I	Number of valid points in reduced grid	Output
C----------------------------------------------------------------------
      INTEGER NX,NY,NH,M,NW,NS,NT,MX,MY
      REAL*4 A(NX,NY),H(NH),PMIN,PMAX,PINC
      REAL*4 JA,IA,X(*),Y(*),Z(*)
*
      REAL R,GRAV(2)
      INTEGER I,J,NF,K,IER
*
      REAL CVI(6),CVO(6)
      COMMON/CONV/CVI,CVO
*
      K = 0
      DO J = 1,MY
         JA = (J-CVO(4))*CVO(6)+CVO(5)
         DO I = 1,MX
            IA = (I-CVO(1))*CVO(3)+CVO(2)
            CALL ZODI102(A,NX,NY,IA,JA,PMIN,PMAX,PINC,
     $      H,NH,NW,NF)
            IF (NF.GT.NT) THEN
               CALL ZODI122(H,NH,R,PMIN,PINC)
               CALL ZODI123(A,NX,NY,IA,JA,NW,PINC,R,GRAV,IER)
               IF (IER.EQ.0) THEN
                  K    = K+1
                  Z(K) = R
                  X(K) = GRAV(1)
                  Y(K) = GRAV(2)
               ENDIF
            ENDIF
         ENDDO
      ENDDO
      WRITE(6,*) K,' intermediate background values found'
      IF (K.GT.0) CALL ZODI124(K,X,Y,Z)
      M = K
      WRITE(6,*) M,' unredondant (valid) background values found'
      DO J=1,M
         WRITE(1,*) X(J),Y(J),1.0,Z(J)
      ENDDO
      END
*<FF>
      SUBROUTINE ZODI124(N,X,Y,Z)
C     ensures no two entries have the same coordinates (arises when
C     background is estimated on areas greater than the mesh width)
      REAL X(*),Y(*),Z(*)
      INTEGER N,I,K,GOOD
*
      GOOD=1
      DO I=2,N
         DO K=1,GOOD
            IF (X(I).EQ.X(K).AND.Y(I).EQ.Y(K)) GOTO 1
         ENDDO
         GOOD=GOOD+1
         X(GOOD)=X(I)
         Y(GOOD)=Y(I)
         Z(GOOD)=Z(I)
1     ENDDO
      N=GOOD
      RETURN
      END
*<FF>
      SUBROUTINE ZODI102 (A,NX,NY,IA,JA,PMIN,PMAX,PINC,H,NH,NW,NF)
C----------------------------------------------------------------------
C Determine an histogram for map values in a disk around the
C specified point.
C
C Arguments:
C	A	R*4	Input data cube				Input
C	NX	I	First dimension of A			Input
C	NY	I	Second dimension of A			Input
C	IA	R*4	First coordinate of specified point	Input
C	JA	R*4	Second coordinate of specified point	Input
C	PMIN	R*4	Minimum accepted value			Input
C	PMAX	R*4	Maximum accepted value			Input
C	PINC	R*4	Histogram step				Input
C	H	R*4	Histogram array				Output
C	NH	I	Size of histogram			Input
C	NW	I	Diameter of smoothing circle		Input
C	NF	I	Number of valid data points		Output
C----------------------------------------------------------------------
      INTEGER NX,NY,NH,NW,NF
      REAL*4 A(NX,NY),H(NH),PMIN,PMAX,PINC,IA,JA
*
      INTEGER I,J,K,IMIN,IMAX,JMIN,JMAX
      REAL R
*
      REAL CVI(6),CVO(6)
      COMMON/CONV/CVI,CVO
*
* Zero Histogram
      DO I = 1,NH
         H(I) = 0.0
      ENDDO
      NF = 0
*
* Build it
      JMIN = MAX(1,INT(JA-0.5*NW)+1)
      JMAX = MIN(NY,INT(JA+0.5*NW)-1)
      IMIN = MAX(1,INT(IA-0.5*NW)+1)
      IMAX = MIN(NX,INT(IA+0.5*NW)-1)
      DO J = JMIN,JMAX
         DO I = IMIN,IMAX
            R = A(I,J)
            IF (R.GE.PMIN .AND. R.LT.PMAX) THEN
               K = INT((R-PMIN)/PINC)+1
               H(K) = H(K)+1.0
               NF = NF+1
            ENDIF
         ENDDO
      ENDDO
      END
*<FF>
      SUBROUTINE ZODI122 (H,NH,R,PMIN,PINC)
C----------------------------------------------------------------------
C Determine maximum in an histogram.
C
C Arguments:
C	H	R*4	Histogram array				Input
C	NH	I	Size of histogram			Input
C	R	R*4	Coordinate of histogram maximum		Output
C	PMIN	R*4	Minimum accepted value			Input
C	PMAX	R*4	Maximum accepted value			Input
C	PINC	R*4	Histogram step				Input
C	NW	I	Radius of smoothing circle		Input
C----------------------------------------------------------------------
      INTEGER NH,I,N
      REAL H(NH),R,PMIN,PINC,T
*
      T = H(1)
      N = 1
      DO I=2,NH
         IF (H(I).GT.T) THEN
            N = I
            T = H(I)
         ENDIF
      ENDDO
      R = PMIN+(N-0.5)*PINC
      END
*<FF>
      SUBROUTINE ZODI123(A,NX,NY,IA,JA,NW,PINC,R,GRAV,IER)
C----------------------------------------------------------------------
C Determine the center of gravity of pixels belonging to the histogram
C bin of value R (width PINC) in map A, inside circle of diameter
C NW centered at IA JA
C
C Arguments:
C	A	R*4	Input data cube				Input
C	NX	I	First dimension of A			Input
C	NY	I	Second dimension of A			Input
C	IA	R*4	First coordinate of specified point	Input
C	JA	R*4	Second coordinate of specified point	Input
C	NW	I	Diameter of smoothing circle		Input
C	PINC	R*4	Histogram step				Input
C	R	R*4	Coordinate of histogram maximum		Output
C	GRAV(2)	R*4	gravity center coordinates		Output
C	IER	I	Logical error flag			Output
C----------------------------------------------------------------------
      INTEGER NX,NY,NW,IER
      REAL*4 A(NX,NY),PINC,R,IA,JA,GRAV(2)
*
      INTEGER I,J,IMIN,IMAX,JMIN,JMAX,NGRAV
      REAL*4 RMIN,RMAX
*
      REAL CVI(6),CVO(6)
      COMMON/CONV/CVI,CVO
* Zero
      GRAV(1)=0
      GRAV(2)=0
      NGRAV=0
* Build it
      JMIN = MAX(1,INT(JA-0.5*NW)+1)
      JMAX = MIN(NY,INT(JA+0.5*NW)-1)
      IMIN = MAX(1,INT(IA-0.5*NW)+1)
      IMAX = MIN(NX,INT(IA+0.5*NW)-1)
      RMIN=R-PINC*0.5
      RMAX=RMIN+PINC
      DO J = JMIN,JMAX
         DO I = IMIN,IMAX
            IF (A(I,J).GE.RMIN .AND. A(I,J).LE.RMAX) THEN
               GRAV(1)=GRAV(1)+FLOAT(I)
               GRAV(2)=GRAV(2)+FLOAT(J)
               NGRAV=NGRAV+1
            ENDIF
         ENDDO
      ENDDO
      IF (NGRAV.GT.0) THEN
         GRAV(1)=GRAV(1)/NGRAV
         GRAV(2)=GRAV(2)/NGRAV
         IER=0
         RETURN
      ENDIF
      IER=1
      END
*<FF>
      SUBROUTINE ZODI013(M,X,Y,Z,W,IW,B,NX,NY,BVAL,IER)
      INTEGER M,NX,NY,IER
      REAL X(*),Y(*),Z(*),W(*),IW(*),B(NX,NY),BVAL
*
      INTEGER NCP
      LOGICAL LOGRAN(7), ERROR
      EXTERNAL DUMF
      DATA LOGRAN /7*.FALSE./, NCP/4/
*
      CALL GRIDINI(NX,1.D0,1.D0,1.D0, NY,1.D0,1.D0,1.D0)
      CALL GRIDSET(NCP,BVAL,LOGRAN)
      CALL GRIDRAN(X,Y,Z,M,W,IW,B,DUMF,DUMF,DUMF,ERROR)
      IF (ERROR) IER = 10
      END
*
      SUBROUTINE DUMF
      END
*<FF>
      SUBROUTINE ZODI004 (A,B,NX,NY,BVAL,EVAL)
      INTEGER NX,NY
      REAL A(NX,NY),B(NX,NY),BVAL,EVAL
      INTEGER I,J
*
      DO J=1,NY
         DO I=1,NX
            IF (ABS(A(I,J)-BVAL).GT.EVAL .AND.
     $      ABS(B(I,J)-BVAL).GT.EVAL) THEN
               B(I,J) = A(I,J) - B(I,J)
            ELSE
               B(I,J) = BVAL
            ENDIF
         ENDDO
      ENDDO
      END
