C++*********************************************************************
C
C FALB
C                   PROMPTS                       JAN 02 ARDEAN LEITH
C                   OPFILEC                       FEB 03 ARDEAN LEITH
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
C  B-spline interpolation introduced 09/21/89
C  Restriction of the interpolation field 10/13/89
C  Subtraction of one image 03/27/91   -  FALB
C  Quadratic interpolation used as an option. 06/24/91
C  Scratch file on the disk   08/01/91
c
c         FALB
c         FALB_P(BUF,ILIST,NSAM,NROW,LSAM,LROW,NIMA,
c         ANG(RKK,MODE)
c         ENER(CIRC,LCIRC,NRING,NUMR,MODE)
c         ALPRBS(NUMR,NRING,LCIRC,MODE)
c         ALRQ
c         UPDTC(CIRC1,CIRC2,LCIRC,NRING,NUMR,TOT,MAXRIN,IS)
c         OUTRNG
c         CROSRNG
c         FOURING(CIRC,LCIRC,NUMR,NRING,E,MODE)
c         LOG2(N)
c         PRB1D(B,NPOINT,POS)
c         FFTR_D(X,NV)   
c         FFTC_D(BR,BI,LN,KS)
C
C IMAGE_PROCESSING_ROUTINE
C
C23456789012345678901234567890123456789012345678901234567890123456789012
C--*********************************************************************

        SUBROUTINE FALB

        INCLUDE 'CMBLOCK.INC'
        INCLUDE 'CMLIMIT.INC'   

        INTEGER, ALLOCATABLE, DIMENSION(:) :: ILIST
        INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NUMR

        COMMON  /F_SPEC/  FINPAT,NLET,FINPIC
        CHARACTER*80  FINPIC, FINPAT

        INTEGER      MAXRIN
        CHARACTER*1  MODE,NULL,ASK

        DATA  INPIC/77/

        NILMAX = NIMAX
        NULL=CHAR(0)

        ALLOCATE (ILIST(NILMAX), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'AP RA, ILIST',IER)
           RETURN
        ENDIF

C       ASK FOR DATA FILE
        CALL FILELIST(.TRUE.,INPIC,FINPAT,NLET,ILIST,NILMAX,NIMA,
     &      'ENTER TEMPLATE FOR 2-D IMAGE SERIES',IRTFLG)
        IF (IRTFLG .NE. 0) RETURN

C       NIMA - TOTAL NUMBER OF IMAGES
        IF (NIMA .GT. 0)  THEN
           WRITE(NOUT,2001) NIMA
2001       FORMAT(' Number of images: ',I5)
        ELSE
           CALL ERRT(100,'NO IMAGES',NDUM)
           DEALLOCATE(ILIST)
           RETURN
        ENDIF
	
C       GET IMAGE SIZE
        CALL  FILGET(FINPAT,FINPIC,NLET,ILIST(1),INTFLG)

        MAXIM = 0
        CALL OPFILEC(0,.FALSE.,FINPIC,INPIC,'O',IFORM,NSAM,NROW,NSLICE,
     &              MAXIM,' ',.FALSE.,IRTFLG)
        CLOSE(INPIC)
        IF (IRTFLG .NE. 0)  THEN
           DEALLOCATE(ILIST)
           RETURN
        ENDIF
 
        CALL  RDPRMI(MR,NR,NOT_USED,'FIRST AND LAST RING')

        IF(MR.LE.0.OR.NR.GE.MIN0(((NSAM-1)/2)*2+1,((NROW-1)/2)*2+1))THEN
           CALL ERRT(31,'OR 2',NE)
           DEALLOCATE(ILIST)
           RETURN
        ENDIF

        CALL  RDPRMI(ISKIP,NDUMP,NOT_USED,'SKIP')
        ISKIP=MAX0(1,ISKIP)
7981    NA=1
        CALL  RDPRMC(ASK,NA,.TRUE.,'(F)ULL OR (H)ALF CIRCLE',NULL,IRT)

        IF (ASK.EQ.'F')  THEN
           MODE='F'
        ELSEIF (ASK.EQ.'H')  THEN
           MODE='H'
        ELSE
           DEALLOCATE(ILIST)
           RETURN
        ENDIF

C       CALL  RDPRMI(JACUP,NDUMP,NOT_USED,
C     &         'Precision of peak location (0..100)')
C       JACUP=MAX0(0,MIN0(100,JACUP))
        JACUP=0

C       FIND TOTAL NUMBER OF RINGS
        NRING=0
        DO    I=MR,NR,ISKIP
           NRING=NRING+1
        ENDDO

        ALLOCATE (NUMR(3,NRING), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'AP RA, NUMR',IER)
           DEALLOCATE (ILIST)
           RETURN
        ENDIF


        NRING=0
        DO    I=MR,NR,ISKIP
           NRING=NRING+1
           NUMR(1,NRING)=I
        ENDDO
 
C       CALCULATION OF ACTUAL DIMENSION OF AN IMAGE TO BE INTERPOLATED
C       2*(NO. OF RINGS)+(0'TH ELEMENT)+2*(MARGIN OF 1)

        NRA=MIN0(((NSAM-1)/2)*2+1,((NROW-1)/2)*2+1,2*NR+3)
        LSAM=NSAM
        LROW=NROW
        NSAM=NRA
        NROW=NRA
        CALL  ALPRBS(NUMR,NRING,LCIRC,MODE)

        MAXRIN=NUMR(3,NRING)
 
        CALL  FALB_P(ILIST,NSAM,NROW,LSAM,LROW,NIMA,
     &               NRING,LCIRC,MAXRIN,JACUP,NUMR,MODE,NOUT)

         WRITE (NOUT,2600)
2600     FORMAT (/ ' ',72('-')//,
     &             ' ','ROTATIONAL ALIGNMENT -- END OF COMPUTATION',//,
     &             ' ',72('-')/)

        DEALLOCATE (ILIST, NUMR)

        END


C++*********************************************************************
C
C FALB_P.F                ROT FIXED & RANDOMIZED JULY 2000 ARDEAN LEITH
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
C IMAGE_PROCESSING_ROUTINE
C
C        1         2         3         4         5         6         7
C23456789012345678901234567890123456789012345678901234567890123456789012
C--*********************************************************************

	SUBROUTINE  FALB_P(ILIST,NSAM,NROW,LSAM,LROW,NIMA,
     &          NRING,LCIRC,MAXRIN,JACUP,NUMR,MODE,NOUT)

C       BUFIN,ROT,CIROLD,CIRNEW,CIRTMP,TEMP,EC AND DIST ARE AUTOMATIC ARRAYS
	
        REAL, ALLOCATABLE, DIMENSION(:,:) :: X,CIRC 

        PARAMETER  (NLIST=5)

        INTEGER  MAXRIN,MAXRI,NUMR(3,NRING)
        DIMENSION  BUFIN(LSAM),DIST(NIMA)
        DIMENSION  CIROLD(LCIRC),CIRNEW(LCIRC),CIRTMP(LCIRC)
        DOUBLE PRECISION TEMP(MAXRIN,2)
        DOUBLE PRECISION ENER,TOTMIN,SOLD,SNEW,EAV,EC(NIMA)
        DIMENSION  ILIST(NIMA),DLIST(NLIST),ROT(NIMA)

        CHARACTER*80  FINPIC,FINPAT

        COMMON  /F_SPEC/  FINPAT,NLET,FINPIC
        COMMON  /MXR/  MAXRI

        LOGICAL*1    CH_ANG
        CHARACTER*1  MODE
 
        DATA  INPIC/77/,NDOC/55/

        MAXRI = MAXRIN

        LQ  = LROW/2+1
        LR1 = (NROW-1)/2
        LR2 = LQ+LR1
        LR1 = LQ-LR1
        LQ  = LSAM/2+1
        LS1 = (NSAM-1)/2
        LS2 = LQ+LS1
        LS1 = LQ-LS1

        ALLOCATE (X(NSAM,NROW), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'AP RA, X',IER)
           RETURN
        ENDIF

        ALLOCATE (CIRC(LCIRC,NIMA), STAT=IRTFLG)
        IF (IRTFLG.NE.0) THEN 
           CALL ERRT(46,'AP RA, CIRC',IER)
           DEALLOCATE (X)
           RETURN
        ENDIF

        DO K1=1,NIMA

           CALL  FILGET(FINPAT,FINPIC,NLET,ILIST(K1),INTFLAG)
           IF (IRTFLG .NE. 0) THEN
              CALL ERRT(18,'AP RA ',NE)
              DEALLOCATE (X,CIRC)     
              RETURN
           ENDIF

           MAXIM = 0
           CALL OPFILEC(0,.FALSE.,FINPIC,INPIC,'O',IFORM,
     &           NSAMT,NROWT,NSL,MAXIM,'DUMMY',.FALSE.,IRTFLG)
          
           IF (IRTFLG .NE. 0)  THEN
              CALL ERRT(4,'AP RA ',NE)
              DEALLOCATE (X,CIRC)
              RETURN
           ENDIF

           DO K2=LR1,LR2
              CALL  REDLIN(INPIC,BUFIN,LSAM,K2)
              DO K3=LS1,LS2
                 X(K3-LS1+1,K2-LR1+1)=BUFIN(K3)
              ENDDO
           ENDDO
       
           CLOSE(INPIC)

           CALL ALRQ(X,NSAM,NROW,NUMR,CIRC(1,K1),LCIRC,NRING,MODE,K1)
           CALL FOURING(CIRC(1,K1),LCIRC,NUMR,NRING,EC(K1),MODE)

        ENDDO

C       BUILD FIRST AVERAGE

C       TWO ESTIMATION OF INITIAL AVERAGE ARE USED
C       ONLY ONE !!!  11/06/91

C          DIST IS USED HERE FOR THE RANDOM CHOOSING OF IMAGES

           DO IMI=1,NIMA
              DIST(IMI) = 0.0
           ENDDO

           CALL RANDOM_NUMBER(CIID)
           IMI = MIN0(NIMA,MAX0(1,INT(CIID*NIMA+0.5)))

           DO I=1,LCIRC
              CIROLD(I) = CIRC(I,IMI)
           ENDDO
           ROT(IMI)  = 1.0
           DIST(IMI) = 1.0

C          write(nout,*) 'rot(',imi,'):',rot(imi),ciid,nima,dist(imi)
           DO KTN=2,NIMA

804           CALL RANDOM_NUMBER(CIID) 
              M = MIN0(NIMA,MAX0(1,INT(CIID*(NIMA-KTN+1)+0.5)))

              IMI = 0
              DO I=1,NIMA
                 IF (DIST(I) .NE. 1.0)  THEN
                    IMI = IMI + 1
                    IF (IMI. EQ. M)  GOTO 810
                 ENDIF
809           CONTINUE
              ENDDO
              GOTO  804

810           IMI       = I
              DIST(IMI) = 1.0

              CALL CROSRNG(CIROLD,CIRC(1,IMI),LCIRC,NRING,TEMP,
     &             TEMP(1,2),MAXRIN,JACUP,NUMR,TOTMIN,TOT,MODE)
              ROT(IMI) = TOT
C             write(nout,*) 'rot(',imi,'):',rot(imi),ciid,m,dist(imi),i
              CALL UPDTC(CIROLD,CIRC(1,IMI),LCIRC,NRING,NUMR,
     &                   TOT,MAXRIN,KTN)
           ENDDO

           SOLD = ENER(CIROLD,LCIRC,NRING,NUMR,MODE)
           WRITE(NOUT,2037)  SOLD
2037       FORMAT(' Random approximation #1,  Squared sum =',1PE12.5)

           DO I=1,LCIRC
              CIRNEW(I) = CIROLD(I)
           ENDDO
                 
C       PRINT  *,SOLD*FLOAT(NIMA)/(NIMA-1)
C       PRINT  2001,(ANG(ROT(J),MODE),J=1,NIMA)
2001    FORMAT(8(1X,F8.3))

C       ITERATIONS TO GET BETTER APPROXIMATION

        ITER=0
901     CONTINUE

        ITER   = ITER+1
        CH_ANG = .FALSE.

        DO IMI=1,NIMA
           CALL OUTRNG(CIROLD,CIRC(1,IMI),CIRTMP,LCIRC,NRING,
     &                 NUMR,ROT(IMI),MAXRIN,NIMA)

           CALL  CROSRNG(CIRTMP,CIRC(1,IMI),LCIRC,NRING,TEMP,
     &               TEMP(1,2),MAXRIN,JACUP,NUMR,TOTMIN,TOT,MODE)

           IF (ROT(IMI) .NE. TOT)  THEN
              CH_ANG   = .TRUE.
              ROT(IMI) = TOT
           ENDIF
          CALL UPDTC(CIRNEW,CIRC(1,IMI),LCIRC,NRING,NUMR,TOT,MAXRIN,IMI)

        ENDDO

        SNEW=ENER(CIRNEW,LCIRC,NRING,NUMR,MODE)
        WRITE(NOUT,2030)  ITER,SNEW
2030    FORMAT('  Iteration #',I4,'  New squared sum =',1PE12.5)
        IF (SNEW.GE.SOLD .AND. CH_ANG)  THEN
           DO  I=1,LCIRC 
              CIROLD(I)=CIRNEW(I)
           ENDDO
           SOLD = SNEW
           GOTO  901
        ENDIF

        WRITE(NOUT,2001)  (ANG(ROT(J),MODE),J=1,NIMA)
        DLIST(5) = 1.0
        DO   IMI=1,NIMA
C          CALCULATE DISTANCES
           CALL  OUTRNG(CIROLD,CIRC(1,IMI),CIRTMP,LCIRC,NRING,
     &           NUMR,ROT(IMI),MAXRIN,NIMA)

           CALL CROSRNG (CIRTMP,CIRC(1,IMI),LCIRC,NRING,TEMP,
     &          TEMP(1,2), MAXRIN,JACUP,NUMR,TOTMIN,TOT,MODE)

           EAV      = ENER(CIRTMP,LCIRC,NRING,NUMR,MODE)
           DLIST(1) = IMI
           DLIST(2) = ILIST(IMI)
           DLIST(3) = ANG(ROT(IMI),MODE)
           DLIST(4) = EAV+EC(IMI)-2.0*TOTMIN
           CALL SAVD(NDOC,DLIST,NLIST,IRTFLG)
        ENDDO
        CALL  SAVDC
        CLOSE(NDOC)

        DEALLOCATE (X,CIRC)
        END