C++*******************************************************************
C
C AVERG3.F
C
C **********************************************************************
C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM.   AUTHOR: J.FRANK  *
C=* Copyright (C) 1985-2005  Health Research Inc.                      *
C=*                                                                    *
C=* HEALTH RESEARCH INCORPORATED (HRI),                                *   
C=* ONE UNIVERSITY PLACE, RENSSELAER, NY 12144-3455.                   *
C=*                                                                    *
C=* Email:  spider@wadsworth.org                                       *
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 the *
C=* 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 GNU  *
C=* General Public License for more details.                           *
C=*                                                                    *
C=* You should have received a copy of the GNU General Public License  *
C=* along with this program; if not, write to the                      *
C=* Free Software Foundation, Inc.,                                    *
C=* 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.      *
C=*                                                                    *
C **********************************************************************
C      AVERG3(LUN1,XAVG,J,K,L,NSAM,NROW,NSLICE,DS,DR,DL,MOVWAY)
C      
C      FOR FURTHER INFORMATION, SEE COMMENTS IN AVERG3.
C
C  XAVG: RECURSIVE LOCAL AVERAGE OF 3D PICTURE VOLUME
C  X( ): DL PICTURE SLICES STORED
C  X[J,K,L]<>X(MOD(L-1,DL)*NROW+J*NSAM+K)
C  J,K,L PICTURE COORD
C  NSAM,NROW,NSLICE: PICTURE SIZE
C  DS,DR,DL: LOCAL BOX  SIZE
C  MOVWAY= 1: READ IN DL SLICES AND COMPUTE 1ST AVERAGE VALUE FROM SCRATCH
C             IF J,K,L IS IN A BORDER AREA COMPUTE XAVG FOR NEAREST
C             ACTIVE POINT
C        = 2: MOVE RT ONE PIXEL
C        = 3: MOVE LEFT ONE PIXEL
C        = 4: MOVE DOWN RIGHT ONE PIXEL
C        = 5: MOVE DOWN LEFT  ONE PIXEL
C        = 6: MOVE UP, READ IN NEXT SLICE
C
C IMAGE_PROCESSING_ROUTINE
C
C        1         2         3         4         5         6         7
C23456789012345678901234567890123456789012345678901234567890123456789012
C--*******************************************************************

        SUBROUTINE AVERG3(LUN1,XAVG,J,K,L,NSAM,NROW,NSLICE,
     &                   DS,DR,DL,MOVWAY)

        INTEGER DS,DR,DL,HDS,HDR,HDL,EDGFST,EDGLST,EDGLFT,EDGRT
        INTEGER EDGDW,EDGUP
        COMMON /UNITS/LUN,NIN,NOUT,NECHO,IFOUND,NPROC,NDAT
        DIMENSION X(1)
        COMMON X
        SAVE

        IF (MOVWAY .EQ. 1) THEN        
C          INITIALIZE PARM AND XAVG
           HDR=DR/2
           HDS=DS/2
           HDL=DL/2
           EDGFST=HDR+1
           EDGLST=NROW-HDR
           EDGLFT=HDS+1
           EDGRT=NSAM-HDS
           EDGDW =HDL+1
           EDGUP=NSLICE-HDL

C          CHECK IF IN BORDER
           XAVG=0.
           JT=J
           KT=K
           KL=L
           IF(J.LT.EDGFST)JT=EDGFST
           IF(J.GT.EDGLST)JT=EDGLST
           IF(K.LT.EDGLFT)KT=EDGLFT
           IF(K.GT.EDGRT)KT=EDGRT
           IF(L.LT.EDGDW) KL=EDGDW
           IF(L.GT.EDGUP) KL=EDGUP
C          READ IN DL SLICES AND COMPUTE AVERAGE FROM SCRATCH
           DO    LL=1,DL
              LSTART=(LL-1)*NSAM*NROW
              DO    JJ=1,NROW
                 NL=MOD(LL-1,DL)*NROW+JJ
                 NI=LSTART+(JJ-1)*NSAM
                 CALL  REDLIN(LUN1,X(NI+1),NSAM,NL)
              ENDDO
           ENDDO

           DO  LL=KL-HDL,KL+HDL
              LSTART=MOD(LL-1,DL)*NSAM*NROW
              DO  JJ=JT-HDR,JT+HDR
                 NI=LSTART+(JJ-1)*NSAM
                 DO  KK=KT-HDS,KT+HDS
                    XAVG=XAVG+X(NI+KK)
                 ENDDO
              ENDDO
           ENDDO
           XAVG=XAVG/FLOAT(DR*DS*DL)
           RETURN

        ELSEIF (MOVWAY .EQ. 2) THEN 
C          MOVE RT ONE PIXEL, UPDATE XAVG
           JT=J
           LT=L
           IF(K.LE.EDGLFT.OR.K.GT.EDGRT)RETURN
           IF(J.LT.EDGFST)JT=EDGFST
           IF(J.GT.EDGLST)JT=EDGLST
           IF(L.LT.EDGDW) LT=EDGDW
           IF(L.GT.EDGUP) LT=EDGUP
           CORECT=0.
           DO    LL=LT-HDL,LT+HDL
              LSTART=MOD(LL-1,DL)*NROW*NSAM
              DO  JJ=JT-HDR,JT+HDR
                 NI=LSTART+(JJ-1)*NSAM
                 CORECT=CORECT+X(NI+K+HDS)-X(NI+K-HDS-1)
              ENDDO
           ENDDO
           XAVG=CORECT/FLOAT(DR*DS*DL)+XAVG
           RETURN

        ELSEIF (MOVWAY .EQ. 3) THEN
C          MOVE LEFT ONE PIXEL UPDATE XAVG
           JT=J
           LT=L
           IF(K.LT.EDGLFT.OR.K.GE.EDGRT)RETURN
           IF(J.LT.EDGFST)JT=EDGFST
           IF(J.GT.EDGLST)JT=EDGLST
           IF(L.LT.EDGDW) LT=EDGDW
           IF(L.GT.EDGUP) LT=EDGUP
           CORECT=0.
           DO    LL=LT-HDL,LT+HDL
              LSTART=MOD(LL-1,DL)*NROW*NSAM
              DO  JJ=JT-HDR,JT+HDR
                 NI=LSTART+(JJ-1)*NSAM
                 CORECT=CORECT+X(NI+K-HDS)-X(NI+K+HDS+1)
              ENDDO
           ENDDO
           XAVG=CORECT/FLOAT(DR*DS*DL)+XAVG
           RETURN

        ELSEIF (MOVWAY .EQ. 4) THEN
C          J INCREMENTED ONE
           IF(J.LE.EDGFST.OR.J.GT.EDGLST)RETURN
           LT=L
           KT=K
           IF(K.LT.EDGLFT)KT=EDGLFT
           IF(K.GT.EDGRT)KT=EDGRT
           IF(L.LT.EDGDW) LT=EDGDW
           IF(L.GT.EDGUP) LT=EDGUP
           CORECT=0.
           DO    LL=LT-HDL,LT+HDL
              LSTART=MOD(LL-1,DL)*NROW*NSAM
              DO  KK=KT-HDS,KT+HDS
                 NI=LSTART+KK+(J-HDR-2)*NSAM
                 NO=LSTART+KK+(J-1+HDR)*NSAM
                 CORECT=CORECT+X(NO)-X(NI)
              ENDDO
           ENDDO
           XAVG=XAVG+CORECT/FLOAT(DS*DR*DL)
           RETURN

        ELSEIF (MOVWAY .EQ. 5) THEN
C          J DECREMENTED ONE
           IF(J.LT.EDGFST.OR.J.GE.EDGLST)RETURN
           LT=L
           KT=K
           IF(K.LT.EDGLFT)KT=EDGLFT
           IF(K.GT.EDGRT)KT=EDGRT
           IF(L.LT.EDGDW) LT=EDGDW
           IF(L.GT.EDGUP) LT=EDGUP
           CORECT=0.
           DO    LL=LT-HDL,LT+HDL
              LSTART=MOD(LL-1,DL)*NROW*NSAM
              DO  KK=KT-HDS,KT+HDS
                 NI=LSTART+KK+(J+HDR)*NSAM
                 NO=LSTART+KK+(J-HDR-1)*NSAM
                 CORECT=CORECT+X(NO)-X(NI)
              ENDDO
           ENDDO
           XAVG=XAVG+CORECT/FLOAT(DS*DR*DL)
           RETURN

        ELSEIF (MOVWAY .EQ. 6) THEN
C          L INCREMENTED ONE
           IF(L.LE.EDGDW.OR.L.GT.EDGUP)  RETURN
           JT=J
           IF(J.LT.EDGFST)JT=EDGFST
           IF(J.GT.EDGLST)JT=EDGLST
           KT=K
           IF(K.LT.EDGLFT)KT=EDGLFT
           IF(K.GT.EDGRT)KT=EDGRT
           CORECT=0.0
           LOLD=L-HDL-1
           LSTART=MOD(LOLD-1,DL)*NSAM*NROW
           DO    JJ=JT-HDR,JT+HDR
              NI=(JJ-1)*NSAM+LSTART
              DO  KK=KT-HDS,KT+HDS
                 CORECT=CORECT-X(NI+KK)
              ENDDO
           ENDDO

           DO    JJ=1,NROW
              NL=L*NROW+JJ
              NI=LSTART+(JJ-1)*NSAM
              CALL  REDLIN(LUN1,X(NI+1),NSAM,NL)
           ENDDO

           DO    JJ=JT-HDR,JT+HDR
              NI=(JJ-1)*NSAM+LSTART
              DO  KK=KT-HDS,KT+HDS
                 CORECT=CORECT+X(NI+KK)
              ENDDO
           ENDDO
           XAVG=XAVG+CORECT/FLOAT(DS*DR*DL)

           RETURN
        ENDIF
        END