C===========================================================================
C Copyright (C) 1995-2005 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===========================================================================
      PROGRAM SMOOTH
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program SMOOTH			version 3.30	850729
C  K. Banse      			ESO - Garching
C
C.KEYWORDS
C  smoothing, average, median, template matching
C
C.PURPOSE
C  (1) filter an image with a Gaussian filter
C  (2) filter an image with a digital filter (given values)
C  (3) filter an image with a digital filter (kernel = another frame)
C  (4) smooth an image by a) averaging in the neighborhood
C                         b) taking the median in the neighbourhood
C  (5) expand an image
C  (6) filter an image with a maximum filter
C  (7) filter an image with a minimum filter	FILACT set to these values...
C
C.ALGORITHM
C  (1) get the radius in x and y + threshold for the averaging/median algorithm
C      from INPUTR and calculate new points via:
C      a) B(n) = SUM A(j), j in neighborhood of given radius around n, if
C         abs(SUM - A(n)) > abs(threshold*A(n)), else B(n) = A(n)
C         or if
C         abs(SUM - A(n)) > threshold , else B(n) = A(n)
C      b) B(n) = median of A(j), j in neighborhood of given radius around n, if
C         abs(median - A(n)) > abs(threshold*A(n)), else B(n) = A(n)
C         or if
C         abs(median - A(n)) > threshold , else B(n) = A(n)
C
C  (2) get the weights for the template from INPUTR or IN_B and calculate new
C      points via:
C      B(n) = SUM A(j)*W(l), j in neighborhood around n, l according to template
C      numbering: 1 ................ nw   (nw = 2*radius + 1)
C  	          nw+1 ........... 2*nw
C	          .....................
C	          (nw-1)*nw+1 ... nw*nw
C
C.INPUT/OUTPUT
C  the following keys are used:
C
C  IN_A/C/1/60          input frame
C  OUT_A/C/1/60         output frame
C  INPUTC/C/1/2         (1:1) = Y, include center pixel in calculation
C                  	      = N, do not include center pixel
C                  	(2:2) = A, use absolute threshold
C                  	      = R, use relative threshold with pixel value
C                  	      = Q, use rel. threshold with average/median value
C                  	      = F, use fixed threshold-interval
C                  	      = Z, use a+b*SQRT(median) as threshold
C  ACTION/C/1/5         (1:1) = A for averaging, M for median,
C				D for digital, G for Gaussian
C				K for maximum filter
C				L for minimum filter
C				E for just expanding the frame
C  			(2:2) = L for looping, N for not looping
C  			(3:3) = D for immediate display, N for no display
C  			(4:4) = E for expansion of input frame
C  			(5:5) = W for auxiliary window, N for not
C			(5:5) not implemented yet!
C  INPUTR/R/1/4         radius ( >0 )in x + y,
C                  	threshold(s) for averaging algorithm
C  INPUTC/C/3/60        area for filtering,
C			if set to +, complete frame is filtered
C                  	[...:...] or CURSOR or table_name
C  INPUTC/C/63/2        :x or :y to indicate that 1-dim digital filter
C  IN_B/C/1/60		name of filter image or template
C  INPUTR/R/1/40	for (2) weights of template
C			for (3) xmean,xsigma,ymean,ysigma
C  INPUTI/I/1/2		x-, y-radius for Gaussian filter image (must be odd)
C  INPUTI/I/1/4		xlef,xri,ylef,yri expansion pixels for EXPANSION
C 
C
C.VERSIONS
C  051014	last modif
C
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER      AREAFL,IDUM(1),IAV,IOMODE
      INTEGER      MINSZX,MINSZY,FF,ZFLAG
      INTEGER      HIOFF,HOOFF,HISIZ,HOSIZ,HY,LINCO,REMAIN
      INTEGER*8    PNTRHI,PNTRHO,VPNTR
      INTEGER      HACK,RETBUF(3)
      INTEGER      N,NAXIS,NCOL,NMAL,NOPNTS,NROW,NWORK
      INTEGER      CHUNKI,CHUNKO,K1,K2
      INTEGER*8    PNTRA,PNTRC,PNTRI,PNTRO,PNTRX,PNTRY,FPNTR
      INTEGER      QPIX(2)
      INTEGER      STAT,STAT1,STAT2,SUBDIM,COOFF
      INTEGER      SUBPIX(2,2),SUBOFA(2),SUBOFB(2),RADIUS(2)
      INTEGER      SUBLO(3),SUBHI(3)
      INTEGER      ICUR1(2),ICUR2(2),LOADFL(2),CUTIDX,LOW(3)
      INTEGER      NPIX(3),TPIX(3),RESPIX(3),OPIX(2),KCOOR(4)
      INTEGER      IMNOA,IMNOC,IMNOI,IMNOO,IMNOX,FIMNO,IMNOEX
      INTEGER      TID,TBCLNM(4),ALL,TABNUL(4)
      INTEGER      FPIX(2),FILACT,KLOOP,DAZO,TRUNAX,HCOUNT
      INTEGER      FNAXIS,FSIZE,WW(4),OFFPIX,RESPLZ,INPLZ,COUNT3
      INTEGER      EXNPIX(3),EXBUF(4),OLDMAP,DEBUFL,EXPAZ,OFFEX
      INTEGER      UNIT(1),NULLO,MADRID(1)
C
      CHARACTER*60   FRAMEA,FRAMEC,FRAME,TABLE,ERMES3
      CHARACTER      CUNIT*64,IDENT*72,CBUF*80,ACTION*5,DEFAUL*2
      CHARACTER      TLABL(4)*16,CAREA*60,XYCH(2)*1,MIDUNI*2
      CHARACTER      TMPFIL*14,CACT*8
C
      DOUBLE PRECISION STEP(3),START(3),EXSTA(3)
      DOUBLE PRECISION FSTART(2),FSTEP(2)
C 
      REAL         INPUTR(40),RMIN,RMAX
      REAL         PCUR1(6),PCUR2(6)
      REAL         RDUM(1),KUTS(4),CUTS(4),RBUF(17)
      REAL         DGFIL1(9),DGFIL2(9),DGFIL3(9),DGFIL4(9)
      REAL         DGFIL5(9),DGFIL6(9),DGFIL7(25),DGFIL8(25)
C
      LOGICAL      SELFLG
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C 
      INCLUDE  'MID_INCLUDE:IDIDEV.INC'
      INCLUDE  'MID_INCLUDE:IDIMEM.INC'
C 
      COMMON   /VMR/ MADRID
C
      DATA    SUBLO    /3*1/,     SUBHI  /3*1/
      DATA    TPIX     /3*1/,     LOW    /3*1/
      DATA    NPIX     /3*1/,     KCOOR  /-1,-1,-1,-1/
      DATA    RESPIX   /3*1/,     KUTS   /4*0.0/
      DATA    IDENT    /' '/,     CUNIT  /' '/
      DATA    KLOOP    /0/,       LOADFL /0,0/
      DATA    DAZO     /0/,       FIMNO  /-1/
      DATA    ACTION   /' '/,     COOFF  /0/
      DATA    EXNPIX   /1,1,1/,   EXSTA  /0.,0.,0./
      DATA    IMNOI    /-1/,      IMNOO  /-1/
      DATA    OLDMAP   /0/,       EXPAZ  /1/
      DATA    START    /3*0.D0/,  STEP   /3*1.D0/
C 
      DATA    ERMES3   /' both cursors are off, set one on ... '/
      DATA    TLABL    /'XSTART ','YSTART ','XEND ','YEND '/
      DATA    XYCH     /'x','y'/
      DATA    CACT     /'NYYY?C0 '/
C 
C  data for filters CONE, LAPLACE, LOW_PASS, POINT, CONE5:
      DATA    DGFIL1 /1.,-6.,1.,-6.,36.,-6.,1.,-6.,1./         !CONE
      DATA    DGFIL2 /0.,-1.,0.,-1.,5.,-1.,0.,-1.,0./          !LAPLACE
      DATA    DGFIL3 /1.,2.,1.,2.,5.,2.,1.,2.,1./              !LOW_PASS
      DATA    DGFIL4 /-1.,-2.,-1.,-2.,13.,-2.,-1.,-2.,-1./     !POINT
      DATA    DGFIL5 /1.,1.,1.,1.,8.,1.,1.,1.,1./              !TENT
      DATA    DGFIL6 /-1.,-2.,-1.,-2.,13.,-2.,-1.,-2.,-1./     !SHARP
      DATA    DGFIL7 /0.,0.,1.,0.,0.,0.,2.,2.,2.,0.,           !CONE5
     +                1.,2.,5.,2.,1.,0.,2.,2.,2.,0.,0.,0.,1.,0.,0./
      DATA    DGFIL8 /1.,2.,3.,2.,1.,2.,4.,6.,4.,2.,           !PYRAMID5
     +                3.,6.,9.,6.,3.,2.,4.,6.,4.,2.,1.,2.,3.,2.,1./
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C  set up MIDAS environment + get current MIDAS unit
      CALL STSPRO('SMOOTH')
      FILACT = 0
      AREAFL = 0
      OFFEX = 0
      INPLZ = 0
      RESPLZ = 0
      CALL STKRDC('MID$SESS',1,11,2,IAV,MIDUNI,UNIT,NULLO,STAT)
      TMPFIL(1:) = 'filt'//MIDUNI//'dum.bdf '
C
C  get method + special options and branch according to filter type
      CALL STKRDC('ACTION',1,1,5,IAV,ACTION,UNIT,NULLO,STAT)
      CALL UPCAS(ACTION,ACTION)
      IF (ACTION(2:2).EQ.'L') KLOOP = 1                 !update loop flag
      IF (ACTION(3:3).EQ.'D') LOADFL(1) = -99           !set display flag
      IF ((ACTION(4:4).EQ.'N') .AND.                    !clear expansion flag
     +    (LOADFL(1).NE.-99)) EXPAZ = 0                 !only if no display...
      CALL STKRDC('MID$SPEC',1,1,5,IAV,CAREA,UNIT,NULLO,STAT)
      IF (CAREA(1:5).EQ.'DEBUG') THEN
         DEBUFL = 1
      ELSE
         DEBUFL = 0
      ENDIF
C 
      FF = 0
      PNTRHI = -1                               !to know if it's first time
      IF (ACTION(1:1).EQ.'G') THEN                    !Gaussian
         CALL STKRDI('INPUTI',1,2,IAV,RADIUS,UNIT,NULLO,STAT)
         CALL STKRDR('INPUTR',1,4,IAV,INPUTR,UNIT,NULLO,STAT)
         CALL STKRDC('IN_B',1,1,60,IAV,FRAME,UNIT,NULLO,STAT)
         ZFLAG = 1                            !flag for routine TEMPLB
      ELSE IF (ACTION(1:1).EQ.'D') THEN               !digital
         CALL STKRDC('IN_B',1,1,60,IAV,FRAME,UNIT,NULLO,STAT)
         ZFLAG = 2  
      ELSE IF ((ACTION(1:1).EQ.'K') .OR. (ACTION(1:1).EQ.'L')) THEN
         CALL STKRDI('INPUTI',1,2,IAV,RADIUS,UNIT,NULLO,STAT)
         FILACT = 6
         IF (ACTION(1:1).EQ.'K') THEN
            ZFLAG = 3                            !maximum filter
         ELSE
            ZFLAG = 4                            !minimum filter
         ENDIF
      ELSE
         CALL STKRDC('INPUTC',1,1,2,IAV,DEFAUL,UNIT,NULLO,STAT)
         CALL UPCAS(DEFAUL,DEFAUL)
         IF ((DEFAUL(2:2).EQ.'F') .OR. (DEFAUL(2:2).EQ.'Z')) THEN
            CALL STKRDR('INPUTR',1,2,IAV,INPUTR,UNIT,NULLO,STAT)
            CALL STKRDC('P3',1,1,60,IAV,CBUF,UNIT,NULLO,STAT)
            N = INDEX(CBUF,',')
            CBUF(N:N) = ' '
            N = INDEX(CBUF,',') + 1               !skip over 2. comma
            CBUF(1:) = CBUF(N:)                   !shift to front
C 
            IF (DEFAUL(2:2).EQ.'Z') THEN
               CALL GENCNV(CBUF,2,2,IDUM,INPUTR(3),STEP,IAV)
               IF (IAV.LT.2) CALL STETER
     +                  (18,'Invalid syntax for threshold factors')
            ELSE
               IF (CBUF(1:1).EQ.'<') THEN
                  FF = 1
                  CALL GENCNV(CBUF(3:),2,1,IDUM,RDUM,STEP,IAV)
                  INPUTR(4) = RDUM(1)
               ELSE
                  CALL GENCNV(CBUF,2,1,IDUM,RDUM,STEP,IAV)
                  IF (IAV.LT.1) CALL STETER
     +                     (18,'Invalid syntax for interval limits')
                  INPUTR(3) = RDUM(1)
                  N = INDEX(CBUF,',') + 1
                  IF (CBUF(N:N) .EQ. '>') THEN
                     FF = 2
                  ELSE
                     CALL GENCNV(CBUF(N:),2,1,IDUM,RDUM,STEP,IAV)
                     INPUTR(4) = RDUM(1)
                  ENDIF
               ENDIF
               IF (IAV.LT.1) 
     +            CALL STETER(18,'Invalid syntax for interval limits')
            ENDIF
         ELSE
            CALL STKRDR('INPUTR',1,3,IAV,INPUTR,UNIT,NULLO,STAT)
         ENDIF
      ENDIF
C 
C
C  get area option, input frame + output frame 
      CAREA(1:) = ' '
      CALL STKRDC('INPUTC',1,3,60,IAV,CAREA,UNIT,NULLO,STAT)
      CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,UNIT,NULLO,STAT)
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,UNIT,NULLO,STAT)
      CALL GENEQF(FRAMEA,FRAMEC,IAV)  !check that in, out_frame are different
      IF (IAV.NE.0)
     +   CALL STETER
     +   (2,'input and result frame have to be different...')
C 
      CALL STFOPN(FRAMEA,D_R4_FORMAT,0,F_IMA_TYPE,IMNOA,STAT)
      CALL STDRDI(IMNOA,'NAXIS',1,1,IAV,NAXIS,UNIT,NULLO,STAT)
      CALL STDRDI(IMNOA,'NPIX',1,NAXIS,IAV,NPIX,UNIT,NULLO,STAT)
      CALL STDRDD(IMNOA,'START',1,NAXIS,IAV,START,UNIT,NULLO,STAT)
      CALL STDRDD(IMNOA,'STEP',1,NAXIS,IAV,STEP,UNIT,NULLO,STAT)
      CALL STDRDC(IMNOA,'IDENT',1,1,72,IAV,IDENT,UNIT,NULLO,STAT)
      IF (NAXIS.GT.3) CALL STETER(1,'max. 3-dim frames supported...')
      N = (NAXIS+1) * 16
      CALL STDRDC(IMNOA,'CUNIT',1,1,N,IAV,CUNIT,UNIT,NULLO,STAT)
C 
      CALL STDRDR(IMNOA,'LHCUTS',1,4,IAV,CUTS,UNIT,NULLO,STAT)
      IF (CUTS(1).GE.CUTS(2)) THEN
         CUTIDX = 3
      ELSE
         CUTIDX = 1
      ENDIF
C 
C  check for < or > in fixed interval option (only average+median)
      IF (FF.EQ.1) THEN
         INPUTR(3) = CUTS(3) - 0.1
      ELSE IF (FF.EQ.2) THEN
         INPUTR(4) = CUTS(4) + 0.1
      ENDIF
C
C  determine "real" NAXIS
      TRUNAX = NAXIS
      IF (NAXIS.EQ.3) THEN
         IF (NPIX(3).EQ.1) THEN
            IF (NPIX(2).EQ.1) THEN
               TRUNAX = 1
            ELSE
               TRUNAX = 2
            ENDIF
         ENDIF
      ELSE IF (NAXIS.EQ.2) THEN
         IF (NPIX(2).EQ.1) TRUNAX = 1
      ENDIF
C 
C  only for complete frames do we support 3-dim frames 
      IF (TRUNAX.EQ.3) THEN
         IF (CAREA(1:1).NE.'+') THEN
            CBUF(1:) = 
     +      'max 2-dim frames supported for chosen area option... '
            CALL STETER(23,CBUF)
         ENDIF
C 
         CBUF(1:) = '3-dim input frame - all planes are processed ... '
         CALL STTPUT(CBUF,STAT)
      ENDIF
C
      COUNT3 = 1                          !counter for planes in 3-dim case
      OFFPIX = 0
      LOADFL(2) = D_R4_FORMAT             !use real data format
C 
C -------------------------------------------------------------------
C 
C  construct Gaussian filter
C 
C -------------------------------------------------------------------
C
500   IF (ACTION(1:1).EQ.'G') THEN 
         IF (NAXIS.EQ.1) RADIUS(2) = 0               !force 1-dim filtering...
         IF ((RADIUS(1).LE.0).OR.(RADIUS(2).LT.0))
     +      CALL STETER(4,'invalid radius chosen...')
         IF (RADIUS(2).GT.0) THEN
            FNAXIS = 2
         ELSE
            FNAXIS = 1
         ENDIF
C
         FPIX(1) = 2*RADIUS(1) + 1                  !get no. of points
         FPIX(2) = 2*RADIUS(2) + 1                  !and no. of lines...
         FSIZE = FPIX(1) * FPIX(2)
         FSTART(1) = 0.D0
         FSTART(2) = 0.D0
         FSTEP(1) = 1.D0
         FSTEP(2) = 1.D0
C 
C  check, if we should save the filter frame
         IF (FRAME(1:1).EQ.'+') THEN                 !No, don't keep...
            FRAME = 'dummy_frame '
            IOMODE = F_X_MODE
         ELSE                                        !Yes, keep filter frame...
            IOMODE = F_O_MODE
         ENDIF
         IF (FIMNO.GT.-1) CALL STFCLO(FIMNO,STAT)           !because of loop...
         CALL STIPUT(FRAME,D_R4_FORMAT,IOMODE,F_IMA_TYPE,
     +               FNAXIS,FPIX,FSTART,FSTEP,
     +               'Gaussian filter','none ',FPNTR,FIMNO,STAT)
C
C  here we really build the filter (and normalize it)
         CALL GAUSS(MADRID(FPNTR),FNAXIS,FPIX,FSTART,FSTEP,INPUTR,
     +              4,RMIN,RMAX)
         CALL NORMIT(MADRID(FPNTR),FSIZE)
         IF (IOMODE.EQ.F_O_MODE) THEN                 !if we keep file, get cuts
            RBUF(3) = CUTS(3)
            RBUF(4) = CUTS(4)
            CALL MYMX(MADRID(FPNTR),FSIZE,RBUF(3))
            RBUF(1) = 0.
            RBUF(2) = 0.
            CALL STDWRR(FIMNO,'LHCUTS',RBUF,1,4,UNIT,STAT)
         ENDIF
         FILACT = 1
C 
C -------------------------------------------------------------------
C 
C  digital filter
C 
C -------------------------------------------------------------------
C 
      ELSE IF (ACTION(1:1).EQ.'D') THEN
         CALL UPCAS(FRAME,CBUF)
         RADIUS(1) = 1                   !default to 3*3 matrix with weights
         RADIUS(2) = 1
         FILACT = 9                      !and fixed filter name
         K1 = 9
C
C  predefined filters
         RDUM(1) = 0.                    !we want to normalize filter kernels
         IF (CBUF(1:4).EQ.'CONE') THEN
            IF (CBUF(5:5).EQ.'5') THEN
               DO 800, N=1,25
                  INPUTR(N) = DGFIL7(N)
                  RDUM(1) = RDUM(1) + INPUTR(N)
800            CONTINUE
               RADIUS(1) = 2
               RADIUS(2) = 2
               K1 = 25
            ELSE
               DO 805, N=1,9
                  INPUTR(N) = DGFIL1(N)
                  RDUM(1) = RDUM(1) + INPUTR(N)
805            CONTINUE
            ENDIF
         ELSE IF (CBUF(1:7).EQ.'LAPLACE') THEN
            DO 810, N=1,9
               INPUTR(N) = DGFIL2(N)
               RDUM(1) = RDUM(1) + INPUTR(N)
810         CONTINUE
         ELSE IF (CBUF(1:8).EQ.'LOW_PASS') THEN
            DO 820, N=1,9
               INPUTR(N) = DGFIL3(N)
               RDUM(1) = RDUM(1) + INPUTR(N)
820         CONTINUE
         ELSE IF (CBUF(1:5).EQ.'POINT') THEN
            DO 830, N=1,9
               INPUTR(N) = DGFIL4(N)
               RDUM(1) = RDUM(1) + INPUTR(N)
830         CONTINUE
         ELSE IF (CBUF(1:4).EQ.'TENT') THEN
            DO 840, N=1,9
               INPUTR(N) = DGFIL5(N)
               RDUM(1) = RDUM(1) + INPUTR(N)
840         CONTINUE
         ELSE IF (CBUF(1:5).EQ.'SHARP') THEN
            DO 850, N=1,9
               INPUTR(N) = DGFIL6(N)
               RDUM(1) = RDUM(1) + INPUTR(N)
850         CONTINUE
         ELSE IF (CBUF(1:5).EQ.'PYRAM') THEN
            DO 860, N=1,25
               INPUTR(N) = DGFIL8(N)
               RDUM(1) = RDUM(1) + INPUTR(N)
860         CONTINUE
            RADIUS(1) = 2
            RADIUS(2) = 2
            K1 = 25
C 
C  any weights file
         ELSE IF (CBUF(1:2).NE.'++') THEN
            IF (FIMNO.GT.-1) CALL STFCLO(FIMNO,STAT)
            CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +                  2,FNAXIS,FPIX,FSTART,FSTEP,
     +                  CBUF,CBUF,FPNTR,FIMNO,STAT)
C
C  make sure, we have odd NPIX(i)
            IF (FPIX(1).EQ.1) THEN               !kernel is a column
               RADIUS(1) = 0
            ELSE
               IF (MOD(FPIX(1),2).EQ.0) THEN
                  FPIX(1) = FPIX(1) - 1
                  WRITE(CBUF,30001) XYCH(1),FPIX(1)
                  CALL STTPUT(CBUF,STAT)
               ENDIF
               RADIUS(1) = FPIX(1)/2
            ENDIF
            IF (FNAXIS.EQ.1) THEN                !kernel is a row
               RADIUS(2) = 0
            ELSE
               IF (MOD(FPIX(2),2).EQ.0) THEN
                  FPIX(2) = FPIX(2) - 1
                  WRITE(CBUF,30001) XYCH(2),FPIX(2)
                  CALL STTPUT(CBUF,STAT)
               ENDIF
               RADIUS(2) = FPIX(2)/2
            ENDIF
            FILACT = 3
C 
C  the weights are entered directly
         ELSE
            CALL STKRDC('P3',1,1,80,IAV,CBUF,UNIT,NULLO,STAT)
            K1 = 1
            DO 880, N=1,80                 !get no. of values entered
               IF (CBUF(N:N).EQ.',') K1 = K1 + 1
               IF (CBUF(N:N).EQ.' ') GOTO 888
880         CONTINUE
888         IF (K1.GT.39) K1 = 39
            IF (MOD(K1,2).EQ.0) THEN
               K1 = K1 - 1
               WRITE(CBUF,30006) K1
               CALL STTPUT(CBUF,STAT)
            ENDIF
C 
            CALL STKRDC('INPUTC',1,63,2,IAV,CBUF,UNIT,NULLO,STAT)
            IF (CBUF(1:1).EQ.':') THEN
               IF ((CBUF(2:2).EQ.'Y').OR.(CBUF(2:2).EQ.'y')) THEN
                  RADIUS(1) = 0
                  RADIUS(2) = K1/2
               ELSE
                  RADIUS(1) = K1/2
                  RADIUS(2) = 0
               ENDIF
            ELSE                      !2-dim square kernel
C 
               IF (K1.EQ.25) THEN
                  RADIUS(1) = 2
                  RADIUS(2) = 2
               ELSE IF (K1.NE.9) THEN
                  IF (K1.LT.23) THEN
                     K1 = 9           !3x3 or 5x5 are only possibilities
                  ELSE
                     K1 = 25
                     RADIUS(1) = 2
                     RADIUS(2) = 2
                  ENDIF
                  WRITE(CBUF,30007) K1
                  CALL STTPUT(CBUF,STAT)
               ENDIF
            ENDIF
C 
            CALL STKRDR('INPUTR',1,K1,IAV,INPUTR,UNIT,NULLO,STAT)
            IF (IAV .LT. K1) THEN
               N = K1 - IAV
               WRITE(CBUF,30005) N
               DO 890, N=IAV+1,K1
                  INPUTR(N) = 0.0
890            CONTINUE
            ENDIF
            FILACT = 2
         ENDIF
C 
C  final check of radii
         IF ( (RADIUS(1).LT.0) .OR. (RADIUS(2).LT.0) .OR.
     +        ((RADIUS(1).EQ.0).AND.(RADIUS(2).EQ.0)) )
     +      CALL STETER(6,'invalid no. of weights...')
C 
C  for predefined filters divide by kernel sum
         IF (FILACT.EQ.9) THEN
            RDUM(1) = 1./RDUM(1)
            DO 895, N=1,K1
               INPUTR(N) = INPUTR(N) * RDUM(1)
895         CONTINUE
            FILACT = 2             !we have numbers
         ENDIF
C
C -------------------------------------------------------------------
C
C  just expand input frame => result frame
C
C -------------------------------------------------------------------
C
      ELSE IF (ACTION(1:1).EQ.'E') THEN
         CALL STKRDI('INPUTI',1,4,IAV,EXBUF,UNIT,NULLO,STAT)
         CALL EXPAX(IMNOA,OFFPIX,FRAMEC,EXBUF,STAT)
         IF (STAT.NE.0) CALL STETER(17,'invalid radius...')
 
         CALL STFOPN(FRAMEC,D_OLD_FORMAT,0,F_IMA_TYPE,IMNOC,STAT)
         CALL DSCUPT(IMNOA,IMNOC,' ',STAT)
         GOTO 9900                               !leave straight away ...
C 
C -------------------------------------------------------------------
C 
C  average (smoothing) or median filter
C 
C -------------------------------------------------------------------
C 
      ELSE IF (FILACT.NE.6) THEN
         FILACT = 4
C
C  test, if center pixel should be included
         IF (DEFAUL(1:1).EQ.'Y') THEN
            ALL = 1
         ELSE
            ALL = 0
         ENDIF
         ACTION(2:2) = DEFAUL(2:2)
C
C  get radii + threshold(s) for the averaging algorithm
         RADIUS(1) = NINT(INPUTR(1))
         IF ( (NAXIS.GT.1) .AND. (NPIX(2).GT.1) ) THEN
            RADIUS(2) = NINT(INPUTR(2))
         ELSE
            RADIUS(2) = 0
         ENDIF
         IF ( (RADIUS(1).LT.0) .OR. (RADIUS(2).LT.0) .OR.
     +        ((RADIUS(1)+RADIUS(2)).EQ.0) )
     +      CALL STETER(5,'invalid x- or y-radius entered...')
C
         NOPNTS = RADIUS(2)*2 + 1
         IF (ACTION(1:1).EQ.'M') THEN
            N = (RADIUS(1)*2 + 1) * NOPNTS 
            NWORK = 2 * N                 !2 buffers of filter size
         ELSE
            IF (RADIUS(2).GT.0) THEN
               NWORK = NOPNTS + NPIX(1)
            ELSE
               NOPNTS = 0
               NWORK = 1025               !for max RADIUS(1) of 512 
            ENDIF
            NWORK = NWORK * 2             !we want double precision
         ENDIF
C 
         IF (KLOOP.GT.1) CALL STFCLO(IMNOX,STAT)
         CALL STFCRE('smoothdumx',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +               NWORK,IMNOX,STAT)
         CALL STFMAP(IMNOX,F_X_MODE,1,NWORK,IAV,PNTRX,STAT)
         IF (ACTION(1:1).EQ.'M') THEN
            VPNTR = PNTRX + N
         ELSE
            PNTRY = PNTRX + (NOPNTS*2)                  !for double precision
         ENDIF
      ENDIF
C
C  ++++++++++++++++++++++++++
C
C  build up the result frame
C
C  ++++++++++++++++++++++++++
C 
      MINSZX = RADIUS(1)*2 + 1       !minimum size of area for given filtersize
      MINSZY = RADIUS(2)*2 + 1
      IF (EXPAZ.EQ.0) THEN              !adapt size of result frame
         RESPIX(1) = NPIX(1) - RADIUS(1) - RADIUS(1)
         RESPIX(2) = NPIX(2) - RADIUS(2) - RADIUS(2)
         EXSTA(1) = START(1) + (RADIUS(1)*STEP(1))
         EXSTA(2) = START(2) + (RADIUS(2)*STEP(2))
      ELSE                              !we do expand...
         RESPIX(1) = NPIX(1) 
         RESPIX(2) = NPIX(2)
         EXSTA(1) = START(1)
         EXSTA(2) = START(2)
      ENDIF
      RESPIX(3) = NPIX(3)
      EXSTA(3) = START(3)
C 
C  if it's complete frame + NPIX(2) > MONITPAR(20) - hack it up
C  but we do it only once in  the beginning...
C 
      IF (KLOOP.LE.1) THEN
         NWORK = RESPIX(3) * RESPIX(2) * RESPIX(1)
         CALL STFCRE(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NWORK,IMNOC,STAT)
         CALL STDWRI(IMNOC,'NAXIS',NAXIS,1,1,UNIT,STAT)
         CALL STDWRI(IMNOC,'NPIX',RESPIX,1,NAXIS,UNIT,STAT)
         CALL STDWRD(IMNOC,'START',EXSTA,1,NAXIS,UNIT,STAT)
         CALL STDWRD(IMNOC,'STEP',STEP,1,NAXIS,UNIT,STAT)
         CALL STDWRC(IMNOC,'IDENT',1,IDENT,1,72,UNIT,STAT)
         N = (NAXIS+1) * 16
         CALL STDWRC(IMNOC,'CUNIT',1,CUNIT,1,N,UNIT,STAT)
         CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNIT,STAT)
C 
         HACK = 0                             !get max. lines for buffers
         CALL HACKUP(RESPIX,D_R4_FORMAT,RETBUF)
         NWORK = RETBUF(1)                    !no. of lines from input frame
C 
         IF (CAREA(1:1).EQ.'+') THEN                   !complete frame
            IF ( (NWORK.LT.RESPIX(2)) .OR. (TRUNAX.EQ.3) ) THEN
               IF (TRUNAX.EQ.3) THEN
                  HACK = 2                     !cube ...
               ELSE
                  HACK = 1
                  N = RETBUF(2)
                  IF (RETBUF(3).GT.0) N = N + 1
                  WRITE(CBUF,20010) N
                  CALL STTPUT(CBUF,STAT)
               ENDIF
               CHUNKI = NWORK*RESPIX(1)         !size of memory to allocate
               CALL STFXMP(CHUNKI,D_R4_FORMAT,PNTRHO,STAT)
            ELSE                               !map complete input frame
               CHUNKI = RESPIX(1) * RESPIX(2)
               CALL STFMAP(IMNOC,F_O_MODE,1,CHUNKI,IAV,PNTRC,STAT)
            ENDIF
C 
C  if not complete area, copy first
         ELSE
            CALL STTDIS
     +      ('first we copy input frame to result frame...',0,STAT)
C 
            IF (EXPAZ.EQ.0) THEN             !copy inner window to result frame
               LOW(1) = RADIUS(1) + 1
               LOW(2) = RADIUS(2) + 1
               CALL COPWI(IMNOA,NPIX,LOW,RESPIX,
     +                    IMNOC,RESPIX,SUBLO,NWORK)
            ELSE                             !straight file copy in chunks
               CHUNKI = NWORK * NPIX(1)
               CALL STFCRE('smoothdumi',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +                     CHUNKI,IMNOI,STAT)
               CALL STFMAP(IMNOI,F_X_MODE,1,CHUNKI,IAV,PNTRI,STAT)
               CALL COPYFX(MADRID(PNTRI),IMNOA,IMNOC,NPIX,NWORK)
               CALL STFCLO(IMNOI,STAT)
            ENDIF
            IMNOI = -1                             !so that we know later on
         ENDIF
C 
         IF (EXPAZ.EQ.0) THEN
            WRITE(CBUF,30003) RESPIX(1),RESPIX(2)
            CALL STTPUT(CBUF,STAT)
            INPLZ = NPIX(1) * NPIX(2)
         ELSE
            INPLZ = 0
         ENDIF
         OFFEX = 0
         RESPLZ = RESPIX(1) * RESPIX(2)
      ENDIF
C 
C  now expand frame to include also 2*RADIUS columns + lines
      EXBUF(1) = RADIUS(1)                  !expand input frame
      EXBUF(2) = RADIUS(1)
      EXBUF(3) = RADIUS(2)
      EXBUF(4) = RADIUS(2)
C 
C  *******************************************************
C  *******************************************************
C
C  expand frame - this is a loop, if 3-dim frames involved
C
C  *******************************************************
C  *******************************************************
C 
1200  HCOUNT = 1
      FRAMEA = ' '
      CALL GENTIM(FRAMEA)                        !get current time
C 
      IF (EXPAZ.EQ.1) THEN
         IF (TRUNAX.EQ.3)THEN
            WRITE(CBUF,20000) COUNT3,FRAMEA(1:40)
         ELSE
            WRITE(CBUF,20001) FRAMEA(1:40)
         ENDIF
         CALL STTPUT(CBUF,STAT)
C 
         CALL EXPAX(IMNOA,OFFPIX,TMPFIL,EXBUF,STAT)
         IF (STAT.NE.0) CALL STETER(17,'invalid radius...')
C 
         IF (CAREA(1:1).EQ.'+') THEN                   !complete frame
            IF (HACK.GE.1) THEN       !if HACK >= 1, use STFOPN for input file
               CALL STFOPN(TMPFIL,D_R4_FORMAT,0,F_IMA_TYPE,
     +                     IMNOEX,STAT)
               CALL STDRDI(IMNOEX,'NPIX',1,2,IAV,EXNPIX,UNIT,NULLO,STAT)
               CALL STDRDD(IMNOEX,'START',1,2,IAV,EXSTA,UNIT,NULLO,STAT)
               IF (PNTRHI.EQ.-1) THEN
                  NWORK = (RETBUF(1) + 2*RADIUS(2)) * EXNPIX(1)
                  CALL STFXMP(NWORK,D_R4_FORMAT,PNTRHI,STAT)
               ENDIF
            ELSE
               CALL STIGET(TMPFIL,D_R4_FORMAT,F_I_MODE,
     +                     F_IMA_TYPE,2,               !we keep NAXIS of FRAMEA
     +                     IAV,EXNPIX,EXSTA,STEP, 
     +                     CBUF,CBUF,PNTRA,IMNOEX,STAT)
            ENDIF
         ELSE
            CALL STFOPN(TMPFIL,D_R4_FORMAT,0,F_IMA_TYPE,IMNOEX,STAT)
            CALL STDRDI(IMNOEX,'NPIX',1,2,IAV,EXNPIX,UNIT,NULLO,STAT)
            CALL STDRDD(IMNOEX,'START',1,2,IAV,EXSTA,UNIT,NULLO,STAT)
         ENDIF
C 
      ELSE
         IF (TRUNAX.EQ.3) THEN
            WRITE(CBUF,22000) COUNT3,FRAMEA(1:40)
            CALL STTDIS(CBUF,0,STAT)
         ENDIF
         IMNOEX = IMNOA
         EXNPIX(1) = NPIX(1)
         EXSTA(1) = START(1)
         EXNPIX(2) = NPIX(2)
         EXSTA(2) = START(2)
C 
         IF (CAREA(1:1).EQ.'+') THEN                   !complete frame
            IF (HACK.GE.1) THEN
               IF (PNTRHI.EQ.-1) THEN
                  NWORK = (RETBUF(1) + 2*RADIUS(2)) * EXNPIX(1)
                  CALL STFXMP(NWORK,D_R4_FORMAT,PNTRHI,STAT)
               ENDIF
            ELSE
               N = EXNPIX(1) * EXNPIX(2)               !map input frame
               CALL STFMAP(IMNOEX,F_I_MODE,1,N,IAV,PNTRA,STAT)
            ENDIF
         ELSE
            CONTINUE                   !nothing to map now for subframes
         ENDIF
      ENDIF
C
C  +++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  get the different area options for all filters
C
C  +++++++++++++++++++++++++++++++++++++++++++++++++++++++
C 
C  if we are here from looping, we already have the area option
      IF (KLOOP.GT.1) THEN
         NMAL = 0
         IF (AREAFL.EQ.2) THEN          !it's the 2. time we're up here...
            COOFF = 0
            GOTO 1500                   !continue with cursor
         ELSE IF (AREAFL.EQ.3) THEN 
            GOTO 2100                   !continue with table
         ELSE
            NMAL = 1
            IF (AREAFL.EQ.3) THEN 
               GOTO 3000                !single subframe (fixed coords)
            ELSE
               OPIX(1) = EXNPIX(1)  
               OPIX(2) = EXNPIX(2)
               GOTO 5000                !full image
            ENDIF
         ENDIF
      ENDIF
C
C  a) complete input image - only here do we support 3-dim files
      IF (CAREA(1:1).EQ.'+') THEN
         AREAFL = 0                     !use complete input frame...
         PNTRO = PNTRC                  !filter image is output image
         PNTRI = PNTRA                  !use full input image
         OPIX(1) = EXNPIX(1)            !dimensions of expanded frame
         OPIX(2) = EXNPIX(2)
         GOTO 5000
      ENDIF
C 
C  b) subimage defined by coordinates
      IF (CAREA(1:1).EQ.'[') THEN
         CALL EXTCOO(IMNOA,CAREA,2,SUBDIM,SUBLO,SUBHI,STAT)
         IF (STAT.NE.0) CALL STETER(4,' ')      !stop, if invalid coordinates...
         SUBPIX(1,1) = SUBLO(1)
         SUBPIX(2,1) = SUBLO(2)
         SUBPIX(1,2) = SUBHI(1)
         SUBPIX(2,2) = SUBHI(2)
         AREAFL = 1
         NMAL = 1
         GOTO 3000
      ENDIF
C
C  c) cursor defined subimage(s)
      CALL UPCAS(CAREA(1:6),CUNIT)
      IF (CUNIT(1:6).NE.'CURSOR') GOTO 2000          !branch to table option
C
C  get main control block for ImageDisplay + attach
      IF (DAZO.EQ.0) THEN
         CALL DTOPEN(1,STAT)
         DAZO = 1
      ENDIF
      NMAL = 0                                 !reset counter
      CBUF(1:) = ' '
      CALL CONCHA(QDSPNO,QOVCH,1,0)            !clear overlay channel
      CALL SETCUR(QDSPNO,2,0,0,KCOOR,STAT)     !set up cursor rectangle
C 
C  input cursor positions and check status
1500  CALL GETCUR(CACT,CBUF,
     +            ICUR1,PCUR1(3),PCUR1(5),RDUM(1),STAT1,
     +            ICUR2,PCUR2(3),PCUR2(5),RDUM(1),STAT2)
C
C  check cursor status
      IF ( (STAT1.EQ.0) .AND. (STAT2.EQ.0) ) THEN
         IF (COOFF.EQ.0) THEN
            COOFF = 1
            CALL STTDIS(ERMES3,0,STAT)
            CBUF(1:) = ' '
            GOTO 1500
         ELSE
            GOTO 8000                          !cursor loop terminated ...
         ENDIF
      ELSE
         NMAL = NMAL + 1                                !update counter
         COOFF = 1
      ENDIF
C
C  and get frame pixel no's
      SUBPIX(1,1) = PCUR1(3)
      SUBPIX(2,1) = PCUR1(4)
      SUBPIX(1,2) = PCUR2(3)
      SUBPIX(2,2) = PCUR2(4)
      AREAFL = 2
      GOTO 3000
C
C  d) table defined subimage(s) - for this option we loop...
2000  CALL CLNTAB(CAREA,TABLE,0)
      CALL TBTOPN(TABLE,F_I_MODE,TID,STAT)
      CALL TBIGET(TID,NCOL,NROW,N,N,N,STAT)           !get no. of rows + columns
      DO 2050, N=1,4
         CALL TBLSER(TID,TLABL(N),TBCLNM(N),STAT)
         IF (TBCLNM(N).LE.0) THEN
            CBUF(1:) = 'column :'//TLABL(N)//'is missing ... '
            CALL STETER(7,CBUF)
         ENDIF
2050  CONTINUE
C 
      AREAFL = 3
      NMAL = 0                                         !init counter
C 
C  if we also want to display the filtered windows, we have to attach 
      IF (LOADFL(1).NE.0) THEN
         IF (DAZO.EQ.0) THEN
            CALL DTOPEN(1,STAT)
            DAZO = 1
         ENDIF
         CALL PIXXCV('INIT',IMNOA,RBUF,STAT)
      ENDIF
C
C  here we loop
2100  IF (NMAL.GE.NROW) THEN
         WRITE(CBUF,30010) NMAL
         CALL STTPUT(CBUF,STAT)
         GOTO 8000
      ENDIF
C
C  get next row of values
      NMAL = NMAL + 1
      CALL TBSGET(TID,NMAL,SELFLG,STAT)
      IF (.NOT.SELFLG) GOTO 2100
C 
      CALL TBRGET(TID,'R*4 ',NMAL,4,TBCLNM,RBUF,TABNUL,STAT)
      PCUR1(3) = (RBUF(1)-START(1))/STEP(1) + 1.
      PCUR1(4) = (RBUF(2)-START(2))/STEP(2) + 1.
      PCUR2(3) = (RBUF(3)-START(1))/STEP(1) + 1.
      PCUR2(4) = (RBUF(4)-START(2))/STEP(2) + 1.
      IF (LOADFL(1).NE.0) THEN            !we have to get related screen pixels
         RBUF(11) = PCUR1(3)
         RBUF(12) = PCUR1(4)
         CALL PIXXCV('_RS',0,RBUF(11),STAT)
         ICUR1(1) = NINT(RBUF(15))
         ICUR1(2) = NINT(RBUF(16))
      ENDIF
C
C  get frame pixel no's
      SUBPIX(1,1) = PCUR1(3)
      SUBPIX(2,1) = PCUR1(4)
      SUBPIX(1,2) = PCUR2(3)
      SUBPIX(2,2) = PCUR2(4)
C
C -------------------------------------------------
C
C  common section for option b), c) and d)
C
C -------------------------------------------------
C
C  get no.  of points involved
3000  TPIX(1) = SUBPIX(1,2) - SUBPIX(1,1) + 1
      TPIX(2) = SUBPIX(2,2) - SUBPIX(2,1) + 1
C 
      IF (TPIX(1).LT.MINSZX) THEN
         WRITE(CBUF,30004) NMAL,XYCH(1)
         CALL STTPUT(CBUF,STAT)
         GOTO 7000
      ENDIF
      IF (TPIX(2).LT.MINSZY) THEN
         WRITE(CBUF,30004) NMAL,XYCH(2)
         CALL STTPUT(CBUF,STAT)
         GOTO 7000
      ENDIF
C
C  the memory allocation scheme wastes space for options c and d), but so what..
C 
      OPIX(1) = TPIX(1) + 2*RADIUS(1)              !size we need as input frame
      OPIX(2) = TPIX(2) + 2*RADIUS(2)
      N = OPIX(1) * OPIX(2)
C 
      IF (N.GT.OLDMAP) THEN
         OLDMAP = N                                !save size
         IF (IMNOI.GT.-1) THEN
            CALL STFCLO(IMNOI,STAT)
            CALL STFCLO(IMNOO,STAT)
         ENDIF
         CALL STFCRE('smoothdumi',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +               N,IMNOI,STAT)
         CALL STFMAP(IMNOI,F_X_MODE,1,N,IAV,PNTRI,STAT)
         N = TPIX(1) * TPIX(2)
         CALL STFCRE('smoothdumo',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +               N,IMNOO,STAT)
         CALL STFMAP(IMNOO,F_X_MODE,1,N,IAV,PNTRO,STAT)
      ENDIF
C 
      IF (EXPAZ.EQ.0) THEN            !get expanded subframe from original input
         SUBOFA(1) = SUBPIX(1,1) - RADIUS(1)
         SUBOFA(2) = SUBPIX(2,1) - RADIUS(2)
         LOW(1) = SUBOFA(1) + OPIX(1) - 1
         LOW(2) = SUBOFA(2) + OPIX(2) - 1
         IF ( (SUBOFA(1).LT.1) .OR. 
     +        (SUBOFA(2).LT.1) .OR.
     +        (LOW(1).GT.NPIX(1)) .OR.
     +        (LOW(2).GT.NPIX(2))) THEN
            WRITE(CBUF,30002) NMAL
            CALL STTPUT(CBUF,STAT)
            GOTO 7000
         ENDIF
      ELSE                            !get expanded subframe from expanded input
         SUBOFA(1) = SUBPIX(1,1)              !x-pixel to start (-radius+radius)
         SUBOFA(2) = SUBPIX(2,1)              !y-pixel to start (-radius+radius)
      ENDIF
      SUBOFB(1) = 1
      SUBOFB(2) = 1
      CHUNKI = EXNPIX(1) * OPIX(2)
      K1 = (SUBOFA(2)-1)*EXNPIX(1) + 1
      CALL STFMAP(IMNOEX,F_I_MODE,K1,CHUNKI,IAV,PNTRA,STAT)
      QPIX(1) = EXNPIX(1)
      QPIX(2) = OPIX(2)
      SUBOFA(2) = 1
      CALL COPYF1(MADRID(PNTRA),QPIX,SUBOFA,OPIX,MADRID(PNTRI),
     +            OPIX,SUBOFB)
      CALL STFUNM(IMNOEX,STAT)
C
5000  IF (AREAFL.GE.2) THEN
         N = 1
         CALL STKWRI('LOG',N,4,1,UNIT,STAT)            !set LOG(4) = 1
      ENDIF
C
C  +++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  and finally do the filtering
C
C  +++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      IF (FILACT.EQ.4) THEN
         IF (HACK.EQ.0) THEN
            PNTRO = PNTRO + OFFPIX
            IF (RADIUS(2).EQ.0) THEN                   !average/median filter
               IF (ACTION(1:1).EQ.'M') THEN
                  CALL AVER1M(ACTION,MADRID(PNTRI),MADRID(PNTRO),
     +                        OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL,
     +                        MADRID(PNTRX),KUTS(3))
               ELSE
                  CALL AVER1A(ACTION,MADRID(PNTRI),MADRID(PNTRO),
     +                        OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL,
     +                        MADRID(PNTRX),KUTS(3))
               ENDIF
            ELSE
               IF (ACTION(1:1).EQ.'M') THEN
                  CALL AVER2M(ACTION,MADRID(PNTRI),MADRID(PNTRO),
     +                        OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL,
     +                        MADRID(PNTRX),MADRID(VPNTR),KUTS(3))
               ELSE
                  IF (ACTION(2:2).EQ.'Q') THEN
                     CALL BVER2A(ACTION,MADRID(PNTRI),MADRID(PNTRO),
     +                           OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL,
     +                           MADRID(PNTRX),MADRID(PNTRY),KUTS(3))
                  ELSE
                     CALL AVER2A(ACTION,MADRID(PNTRI),MADRID(PNTRO),
     +                           OPIX(1),OPIX(2),RADIUS,INPUTR(3),ALL,
     +                           MADRID(PNTRX),MADRID(PNTRY),KUTS(3))
                  ENDIF
               ENDIF
            ENDIF
         ELSE
            HIOFF = 1 + OFFEX
            HOOFF = 1 + OFFPIX
            LINCO = 0
            REMAIN = RETBUF(1)            !we know, that NPIX(2) >= RETBUF(1)
C 
6000        HY = REMAIN + 2*RADIUS(2)
            HISIZ = HY * EXNPIX(1)
            HOSIZ = REMAIN * RESPIX(1)
            CALL STFGET(IMNOEX,HIOFF,HISIZ,IAV,MADRID(PNTRHI),STAT)
            IF (HACK.EQ.1) THEN
               FRAMEA = ' '
               CALL GENTIM(FRAMEA)                        !get current time
               WRITE(CBUF,20020) HCOUNT,FRAMEA(1:40)
               CALL STTDIS(CBUF,0,STAT)
            ENDIF
C 
            IF (RADIUS(2).EQ.0) THEN                   !average/median filter
              IF (ACTION(1:1).EQ.'M') THEN
                  CALL AVER1M(ACTION,MADRID(PNTRHI),MADRID(PNTRHO),
     +                        OPIX(1),HY,RADIUS,INPUTR(3),ALL,
     +                        MADRID(PNTRX),KUTS(3))
               ELSE
                  CALL AVER1A(ACTION,MADRID(PNTRHI),MADRID(PNTRHO),
     +                        OPIX(1),HY,RADIUS,INPUTR(3),ALL,
     +                        MADRID(PNTRX),KUTS(3))
               ENDIF
            ELSE
               IF (ACTION(1:1).EQ.'M') THEN
                  CALL AVER2M(ACTION,MADRID(PNTRHI),MADRID(PNTRHO),
     +                        OPIX(1),HY,RADIUS,INPUTR(3),ALL,
     +                        MADRID(PNTRX),MADRID(VPNTR),KUTS(3))
               ELSE
                  IF (ACTION(2:2).EQ.'Q') THEN
                     CALL BVER2A(ACTION,MADRID(PNTRHI),MADRID(PNTRHO),
     +                           OPIX(1),HY,RADIUS,INPUTR(3),ALL,
     +                           MADRID(PNTRX),MADRID(PNTRY),KUTS(3))
                  ELSE
                     CALL AVER2A(ACTION,MADRID(PNTRHI),MADRID(PNTRHO),
     +                           OPIX(1),HY,RADIUS,INPUTR(3),ALL,
     +                           MADRID(PNTRX),MADRID(PNTRY),KUTS(3))
                  ENDIF
               ENDIF
            ENDIF
C 
            CALL STFPUT(IMNOC,HOOFF,HOSIZ,MADRID(PNTRHO),STAT)
            LINCO = LINCO + REMAIN               !increase line count + check
            REMAIN = NPIX(2) - LINCO
            IF (REMAIN.GT.0) THEN
               IF (REMAIN.GT.RETBUF(1))
     +             REMAIN = RETBUF(1)        !hack up into RETBUF(1)-chunks
               HIOFF = HIOFF + HISIZ - (2*RADIUS(2)*EXNPIX(1))
               HOOFF = HOOFF + HOSIZ
               HCOUNT = HCOUNT + 1
               GOTO 6000                         !loop again
            ENDIF
         ENDIF
C 
      ELSE
         IF (HACK.EQ.0) THEN
            PNTRO = PNTRO + OFFPIX
            IF ((FILACT.EQ.1) .OR. (FILACT.EQ.3)) THEN   
               CALL TEMPLB(ZFLAG,MADRID(PNTRI),MADRID(PNTRO),OPIX,
     +                     MADRID(FPNTR),RADIUS,KUTS(3),STAT)
            ELSE 
               CALL TEMPLB(ZFLAG,MADRID(PNTRI),MADRID(PNTRO),OPIX,
     +                     INPUTR,RADIUS,KUTS(3),STAT)
            ENDIF
            IF (STAT.NE.0)
     +         CALL STETER(5,'x- or y-radius of filter too large...')
C 
         ELSE
            HIOFF = 1 + OFFEX
            HOOFF = 1 + OFFPIX
            LINCO = 0
            REMAIN = RETBUF(1)             !we know, that NPIX(2) >= RETBUF(1)
C 
6600        HY = REMAIN + 2*RADIUS(2)
            HISIZ = HY * EXNPIX(1)
            HOSIZ = REMAIN * RESPIX(1)
            OPIX(2) = HY
            CALL STFGET(IMNOEX,HIOFF,HISIZ,IAV,MADRID(PNTRHI),STAT)
            IF (HACK.EQ.1) THEN
               FRAMEA = ' '
               CALL GENTIM(FRAMEA)                        !get current time
               WRITE(CBUF,20020) HCOUNT,FRAMEA(1:40)
               CALL STTDIS(CBUF,0,STAT)
            ENDIF
C 
            IF ((FILACT.EQ.1) .OR. (FILACT.EQ.3)) THEN   
               CALL TEMPLB(ZFLAG,MADRID(PNTRHI),MADRID(PNTRHO),OPIX,
     +                     MADRID(FPNTR),RADIUS,KUTS(3),STAT)
            ELSE 
               CALL TEMPLB(ZFLAG,MADRID(PNTRHI),MADRID(PNTRHO),OPIX,
     +                     INPUTR,RADIUS,KUTS(3),STAT)
            ENDIF
            IF (STAT.NE.0)
     +         CALL STETER(5,'x- or y-radius of filter too large...')
C 
            CALL STFPUT(IMNOC,HOOFF,HOSIZ,MADRID(PNTRHO),STAT)
            LINCO = LINCO + REMAIN               !increase line count + check
            REMAIN = NPIX(2) - LINCO
            IF (REMAIN.GT.0) THEN
               IF (REMAIN.GT.RETBUF(1)) 
     +            REMAIN = RETBUF(1)            !hack up into RETBUF(1)-chunks
               HIOFF = HIOFF + HISIZ - (2*RADIUS(2)*EXNPIX(1))
               HOOFF = HOOFF + HOSIZ
               HCOUNT = HCOUNT + 1
               GOTO 6600                         !loop again
            ENDIF
         ENDIF
C 
      ENDIF
C
C  if area option used, merge filtered frame back into output frame
      IF (AREAFL.GT.0) THEN
         SUBOFA(1) = 1
         SUBOFA(2) = 1
         IF (EXPAZ.EQ.0) THEN
            SUBOFB(1) = SUBPIX(1,1) - RADIUS(1)
            SUBOFB(2) = SUBPIX(2,1) - RADIUS(2)
         ELSE
            SUBOFB(1) = SUBPIX(1,1)
            SUBOFB(2) = SUBPIX(2,1)
         ENDIF
         CHUNKO = RESPIX(1) * TPIX(2)
         QPIX(1) = RESPIX(1)
         QPIX(2) = TPIX(2)
         K2 = (SUBOFB(2)-1)*RESPIX(1) + 1
         CALL STFMAP(IMNOC,F_IO_MODE,K2,CHUNKO,IAV,PNTRC,STAT)
         SUBOFB(2) = 1
         CALL COPYF1(MADRID(PNTRO),TPIX,SUBOFA,TPIX,MADRID(PNTRC),
     +               QPIX,SUBOFB)
         CALL STFPUT(IMNOC,K2,CHUNKO,MADRID(PNTRC),STAT)
         IF (LOADFL(1).NE.0) THEN       !for display option load subwindows
            LOW(1) = PCUR1(3)
            LOW(2) = PCUR1(4)
            WW(1) = ICUR1(1)
            WW(2) = ICUR1(2)
            WW(3) = SCALX
            WW(4) = SCALY
            CALL LOADWN(LOADFL,IMNOC,RESPIX,LOW,TPIX,WW,CUTS(CUTIDX))
         ENDIF
         CALL STFUNM(IMNOC,STAT)
      ENDIF
C
C  for cursor or table input loop
7000  IF (AREAFL.GE.2) THEN
         N = 0
         CALL STKWRI('LOG',N,4,1,UNIT,STAT)      !reset LOG(4) = 0
         GOTO (1500,2100),AREAFL-1
C 
C  also loop for 3-dim frames (with complete frame option)
      ELSE IF (TRUNAX.EQ.3) THEN
         COUNT3 = COUNT3 + 1
         IF (COUNT3.LE.NPIX(3)) THEN
            OFFPIX = OFFPIX + RESPLZ
            OFFEX = OFFEX + INPLZ
            IF (EXPAZ.NE.0) CALL STFCLO(IMNOEX,STAT)
            GOTO 1200
         ENDIF
      ENDIF
C
C  update descriptors + add history
8000  IF (KLOOP.EQ.0) THEN
         CALL DSCUPT(IMNOA,IMNOC,' ',STAT)
         IF (AREAfl.EQ.0)
     +      CALL STDWRR(IMNOC,'LHCUTS',KUTS,1,4,UNIT,STAT)
         GOTO 9000
      ENDIF
C 
C +++++++++++++++++++++++++++++++++
C 
C  here we loop inside the module
C 
C +++++++++++++++++++++++++++++++++
C 
      IF (KLOOP.EQ.1) THEN
         KLOOP = 99
         IF (ACTION(1:1).EQ.'D') THEN
            WRITE(*,10000)
         ELSE
            WRITE(*,10001)
         ENDIF
      ENDIF
C 
8800  CBUF(1:) = ' '
      CALL STTDIS('> ',99,STAT)
      READ(*,30000) CBUF
      IF ((CBUF(1:1).EQ.'X').OR.(CBUF(1:1).EQ.'x')) THEN
         KLOOP = 0
         GOTO 8000
      ENDIF
      N = INDEX(CBUF,' ')
C 
      IF (ACTION(1:1).EQ.'G') THEN
         IF (CBUF(1:1).NE.'?') THEN
            CALL GENCNV(CBUF(1:N),1,2,RADIUS,RDUM,STEP,IAV)
            IF (IAV.LT.1) THEN
               WRITE(*,10011)
               GOTO 8800
            ELSE IF (IAV.LT.2) THEN
               RADIUS(2) = RADIUS(1)
            ENDIF
         ENDIF
         N = N + 1                       !move to next parameter
         INPUTR(5) = 1.0
         IF ((CBUF(N:N).NE.'?') .AND.
     +       (CBUF(1:1).NE.' ')) THEN
            CALL GENCNV(CBUF(1:N),2,4,IDUM,INPUTR,STEP,IAV)
            IF (IAV.LT.2) THEN
               INPUTR(5) = -1.0
            ELSE IF (IAV.LT.4) THEN
               INPUTR(3) = INPUTR(1)
               INPUTR(4) = INPUTR(2)
            ENDIF
         ENDIF
C 
      ELSE IF (ACTION(1:1).EQ.'D') THEN
         CALL GENCNV(CBUF,2,9,IDUM,INPUTR,STEP,IAV)
         IF (IAV.LT.9) THEN
            FRAME(1:) = CBUF(1:)               !it's a frame
         ELSE
            FRAME(1:2) = '++'
         ENDIF
C 
      ELSE
         IF ((CBUF(N:N).NE.'?') .AND.
     +       (CBUF(1:1).NE.' ')) THEN
            CALL GENCNV(CBUF(1:N),2,3,IDUM,INPUTR,STEP,IAV)
            IF (IAV.LT.3) THEN
               WRITE(*,10010)
               GOTO 8800
            ENDIF
         ENDIF
         N = N + 1                       !move to next parameter
         IF (CBUF(N:N).NE.'?') 
     +      CALL UPCAS(CBUF(N:N+1),DEFAUL)
      ENDIF 
      GOTO 500
C
C  free data
9000  IF (DAZO.NE.0) CALL DTCLOS(QDSPNO)
      CALL STFCLO(IMNOC,STAT)                 !close result frame
      CALL STFCLO(IMNOA,STAT)                 !close input frame
      CALL STFCLO(IMNOEX,STAT)
      IF (DEBUFL.EQ.0) THEN
         IF (EXPAZ.EQ.1) CALL STFDEL(TMPFIL,STAT)
      ENDIF
9900  CALL STSEPI
C
10000 FORMAT(' Enter parameter P3 like in the command syntax,',
     +       ' type X to exit')
10001 FORMAT(' Enter parameters P3 and P4 like in the command syntax,',
     +       ' type X to exit')
10010 FORMAT(' Invalid input, we need 3 real numbers, try again...')
10011 FORMAT(' Invalid input, we need 2 integer numbers, try again...')
20000 FORMAT('expanding plane ',I3,' - ',A)
20001 FORMAT('expanding frame - ',A)
20010 FORMAT('we process frame in ',I2,' chunks ...')
20020 FORMAT('working on chunk ',I2,' - ',A)
22000 FORMAT('plane ',I3,' - ',A)
30000 FORMAT(A)
30001 FORMAT('Warning: No. of ',A,'-weights must be odd => changed to ',
     +       I2,' ...')
30002 FORMAT('invalid subwindow no.',I4,' - so we skip...')
30003 FORMAT('no expansion - NPIX(1,2) of result frame = ',2I6)
30004 FORMAT
     + ('subimage no. ',I4,' too small in ',A,' for given filter size')
30005 FORMAT('missing kernel values (',I2,') are set to 0.0')
30006 FORMAT('Warning: No. of weights must be odd => changed to ',
     +       I2,' ...')
30007 FORMAT
     + ('Warning: No. of weights must be 3x3 or 5x5 => changed to ',
     +       I2,' ...')
30010 FORMAT(I4,' table entries processed ...')
C 
      END

      SUBROUTINE COPWI(IMNOA,NPIXA,BEGA,NPIXW,IMNOB,NPIXB,BEGB,STRIP)
C 
      IMPLICIT NONE
C
      INTEGER   IMNOA,NPIXA(2),BEGA(2),NPIXW(2)
      INTEGER   IMNOB,NPIXB(2),BEGB(2),STRIP
      INTEGER*8 PNTRI,PNTRO
      INTEGER   IMNOI,IMNOO
      INTEGER   CHUNKI,CHUNKO,LOOPLM,RMAIND,STAT,IAV
      INTEGER   IPIX(2),OPIX(2),WPIX(2),ILOW(2),OLOW(2)
      INTEGER   K1,K2,N
      INTEGER   MADRID(1)
C     
      COMMON   /VMR/ MADRID
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C 
      LOOPLM = NPIXW(2) /  STRIP
      RMAIND = NPIXW(2) - (LOOPLM*STRIP)
C 
      CHUNKI = STRIP * NPIXA(1)
      CALL STFCRE('copywork',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            CHUNKI,IMNOI,STAT)
      CALL STFMAP(IMNOI,F_X_MODE,1,CHUNKI,IAV,PNTRI,STAT)
C
      CHUNKO = STRIP * NPIXB(1)
      CALL STFCRE('copzwork',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            CHUNKO,IMNOO,STAT)
      CALL STFMAP(IMNOO,F_X_MODE,1,CHUNKO,IAV,PNTRO,STAT)
C 
      IPIX(1) = NPIXA(1) 
      IPIX(2) = STRIP
      OPIX(1) = NPIXB(1) 
      OPIX(2) = STRIP
      ILOW(1) = BEGA(1)
      ILOW(2) = 1
      OLOW(1) = BEGB(1)
      OLOW(2) = 1
      WPIX(1) = NPIXW(1)
      WPIX(2) = STRIP
C 
      K1 = (BEGA(2)-1)*NPIXA(1) + 1                !full line offsets
      K2 = (BEGB(2)-1)*NPIXB(1) + 1                !full line offsets
      DO 1000, N=1,LOOPLM
         CALL STFGET(IMNOA,K1,CHUNKI,IAV,MADRID(PNTRI),STAT)  !data in
         CALL COPYF1(MADRID(PNTRI),IPIX,ILOW,WPIX,
     +               MADRID(PNTRO),OPIX,OLOW)
         CALL STFPUT(IMNOB,K2,CHUNKO,MADRID(PNTRO),STAT)      !data out
         K1 = K1 + CHUNKI
         K2 = K2 + CHUNKO
1000  CONTINUE
C 
      IF (RMAIND.NE.0) THEN
         CHUNKI = RMAIND * NPIXA(1)
         CHUNKO = RMAIND * NPIXB(1)
         WPIX(2) = RMAIND
         CALL STFGET(IMNOA,K1,CHUNKI,IAV,MADRID(PNTRI),STAT)  !data in
         CALL COPYF1(MADRID(PNTRI),IPIX,ILOW,WPIX,
     +               MADRID(PNTRO),OPIX,OLOW)
         CALL STFPUT(IMNOB,K2,CHUNKO,MADRID(PNTRO),STAT)      !data out
      ENDIF
C
      CALL STFCLO(IMNOI,STAT)
      CALL STFCLO(IMNOO,STAT)
C 
      RETURN
      END
