!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!.COPYRIGHT (C) 1990 European Southern Observatory
!.IDENT smooth.prg
!.AUTHOR Klaus Banse, ESO - Garching
!.KEYWORDS Smoothing, Filtering
!.PURPOSE smooth an image bu box-car filtering
!.VERSION 1.0 Creation 900129
!
!-------------------------------------------------------
!
DEFINE/PARAM P1 ? IMA "Enter input frame: "
DEFINE/PARAM P2 ? IMA "Enter output frame: "
DEFINE/PARAM P3 1.,0. N "Enter radius, threshold: "
DEFINE/PARAM P4 Y C "Enter center_flag, Y(es) or N(o): "
WRITE/KEYWORD IN_A {P1}
WRITE/KEYWORD OUT_A {P2}
WRITE/KEYWORD INPUTR/R/1/2 {P3}
WRITE/KEYWORD HISTORY Smooth
WRITE/KEYWORD DEFAULT {P4}
RUN smooth.exe
To run the application execute the procedure via ![[*]](/images/l2h/icons.gif/cross_ref_motif.gif){kind=icon}
for details).
C @(#)smooth.for 6.3 (ESO-IPG) 1/17/94 17:53:36
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C program SMOOTH version 5.00 900129
C K. Banse ESO - Garching
C
C.KEYWORDS
C smoothing, average
C
C.PURPOSE
C smooth an image by averaging in the neighborhood
C
C.ALGORITHM
C get the names of input + output frames from Keywords IN_A + IN_B,
C get the radius + threshold for the averaging algorithm from the
C Keyword INPUTI and calculate new points via:
C
C B(n) = SUM A(j), j in neighbourhood of given radius around n, if
C SUM - A(n) > threshold*A(n), else B(n) = A(n)
C
C.INPUT/OUTPUT
C the following Keywords are used:
C
C IN_A/C/1/60 input frame
C OUT_A/C/1/60 output frame
C DEFAULT/C/1/1 default = Y, include center pixel in calculation
C of SUM
C INPUTR/R/1/2 radius ( >0 ), threshold ( in [0,1] ) for
C averaging algorithm
C
C-------------------------------------------------------------------
C
PROGRAM SMOOTH
C
IMPLICIT NONE
C
INTEGER NAXISA,NPIXA(2),IAV,STAT,IRAD
INTEGER IMNOA,IMNOC,PNTRA,PNTRC
INTEGER KNULL,KUNIT(1) !KNULL, KUNIT currently only
INTEGER MADRID(1) !placeholders...!
C
CHARACTER*60 FRAMEA,FRAMEC
CHARACTER CUNITA*64,IDENTA*72,DEFAUL*1
C
DOUBLE PRECISION STEPA(2),STARTA(2) !descr. START + STEP are
C !now double precision
C
REAL THRESH,INPUTR(2),CUTS(4)
C
C the following INCLUDE statements may look strange to a Unix user,
C but the ESO preprocessor takes care of the syntax...
C
INCLUDE 'MID_INCLUDE:ST_DEF.INC' !to get the definitions
C !for the ST interfaces...
C
C very important COMMON block VMR for the pointer emulation in Fortran 77 !!!
COMMON /VMR/ MADRID
C
INCLUDE 'MID_INCLUDE:ST_DAT.INC' !to get the data values
C !for the definitions above...
C
DATA IDENTA /' '/ !initialization of IDENTA, CUNITA
DATA CUNITA /' '/ !is necessary for STIGET ...
C
C set up MIDAS environment + enable automatic error abort
C
CALL STSPRO('SMOOTH')
C
C get name of input frame and map it as an image with pixels
C in floating point format
C use the definitions from the include file ST_DEF.INC (e.g. D_R4_FORMAT)
C
CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUNIT,KNULL,STAT)
CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
+ 2,NAXISA,NPIXA,STARTA,STEPA,
+ IDENTA,CUNITA,PNTRA,IMNOA,STAT)
C
C a call to STETER will abort the program ...
C
IF (NAXISA.NE.2)
+ CALL STETER(1,'input frame must have 2 dimensions...')
C
C get name of output frame and create and map it
C Note, that this does not fill in any data values
C
CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,KUNIT,KNULL,STAT)
CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
+ NAXISA,NPIXA,STARTA,STEPA,
+ IDENTA,CUNITA,PNTRC,IMNOC,STAT)
C
C since we use STIGET and STIPUT we cannot have calls of STFGET or STFPUT...
C
C get flag, if center pixel should be included
C get the radius + threshold for the averaging algorithm
C
CALL STKRDC('DEFAULT',1,1,1,IAV,DEFAUL,KUNIT,KNULL,STAT)
CALL STKRDR('INPUTR',1,2,IAV,INPUTR,KUNIT,KNULL,STAT)
IRAD = NINT(INPUTR(1))
THRESH = INPUTR(2)
C
C Now call the subroutine to do the averaging.
C Note, that the address of the input and output array are passed to the
C subroutine using the array MADRID
C
CALL AVERAG(MADRID(PNTRA),MADRID(PNTRC),NPIXA(1),NPIXA(2),
+ IRAD,THRESH,DEFAUL)
C
C update descr. HISTORY with contents of keyword HISTORY + P1, P2, ...
C and copy all non-standard descriptors from input to output frame
C
IDENTA(1:) = ' '
CALL DSCUPT(IMNOA,IMNOC,IDENTA(1:1),STAT)
C
C release files, update keywords and exit
C STSEPI closes all open images, only then are the data values written
C into the result image file
C
CALL STSEPI
C
END
SUBROUTINE AVERAG(A,B,NX,NY,IRAD,THRESH,ALL)
C
C smooth 2-dim image by taking the average of surrounding points,
C if it differs more than THRESH from original value
C
IMPLICIT NONE
C
INTEGER NX,NY,IRAD
INTEGER N,NN,NL,NRL,NP,OFF,OFFA,DIFX,INDX
C
REAL A(1),B(1),THRESH,TEST,COUNT,W
C
CHARACTER ALL*1
C
C initialize COUNT
C to number of pixels in neighbourhood of center
COUNT = FLOAT(IRAD*(IRAD+1)*4) + 1.
C
C first IRAD lines are not treated
C
DO 1000 N=1,IRAD*NX
B(N) = A(N)
1000 CONTINUE
C
C main loop over inner lines
C
DO 5000 NL=IRAD+1,NY-IRAD
OFF = (NL-1)*NX
C
C leave first IRAD values of each line as they are
C
DO 2000 N=OFF+1,OFF+IRAD
B(N) = A(N)
2000 CONTINUE
C
C average in the inner part
C
DO 4000 NP=IRAD+1,NX-IRAD
INDX = OFF + NP !pixel to work on
W = A(INDX) !save original pixel value
TEST = 0.0 !clear sum
C !sum over IRAD lines above
DO 3000 NRL=1,IRAD !and below current line
DIFX = NRL * NX
OFFA = INDX - DIFX - IRAD
DO 2500 NN=OFFA,OFFA+IRAD*2
TEST = TEST + A(NN)
2500 CONTINUE
C
OFFA = OFFA + (2*DIFX)
DO 2800 NN=OFFA,OFFA+IRAD*2
TEST = TEST + A(NN)
2800 CONTINUE
3000 CONTINUE
C
DO 3100 NN=INDX-IRAD,INDX+IRAD !sum over current line
TEST = TEST + A(NN)
3100 CONTINUE
C
C compare average of neighbourhood with pixel in question
C
IF (ALL.NE.'Y') THEN
TEST = (TEST-W) / (COUNT-1.)
ELSE
TEST = TEST/COUNT
ENDIF
C
IF (ABS (TEST-W).GT.ABS (THRESH*W)) THEN
B(INDX) = TEST
ELSE
B(INDX) = W
ENDIF
4000 CONTINUE
C
C leave last IRAD values of each line as they are
C
DO 4100 NN=OFF+NX-IRAD+1,OFF+NX
B(NN) = A(NN)
4100 CONTINUE
C
5000 CONTINUE
C
C end of main loop - last lines are not treated
C
DO 6000 N=(NY-IRAD)*NX+1,NX*NY
B(N) = A(N)
6000 CONTINUE
C
C that's it folks
C
RETURN
END