next up previous contents index
Next: An Application in C Up: Examples Using Image Data Previous: Examples Using Image Data

An Application in FORTRAN

The following source of the program SMOOTH serves as an example of an application program   ready to be run inside the MIDAS environment. This program smooths an image with a box filter.
This program is tied into MIDAS via the following procedure, smooth.prg:
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!.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
Midas 001> @@ smooth inframe outframe rad,thresh centr_flag

The program is written in standard Fortran 77 code with the exceptions which are taken care of by the ESO provided preprocessor. This preprocessor is not needed on a VAX/VMS machine (see chapter [*] 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


Send comments to web@eso.org
Last update: 1998-10-23