C ++******************************************************************** C C HALI_P.F C OPFILEC 02/24/03 al C FINPAT PARAMETER 06/18/08 al C C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2008 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 HALI_P C C PURPOSE: C C PARAMETERS: C C23456789012345678901234567890123456789012345678901234567890123456789012 C*********************************************************************** SUBROUTINE HALI_P(ILIST,NSAM,NROW,LSAM,LROW,NIMA,NRING, & LCIRC,MAXRIN,JACUP,NUMR,NKMAX,MAXIT,MODE,MIRROR,NORM,NOUT, & FINPAT,NLET) C BUFIN,ROT,CIRNEW,TEMP,DIST,EC,CIROLD,CIRSEED,ES,E,IP,IQ C ARE AUTOMATIC ARRAYS INCLUDE 'CMLIMIT.INC' REAL, ALLOCATABLE, DIMENSION(:,:) :: X,CIRC PARAMETER (NLIST=5) INTEGER :: MAXRIN,MAXRI,NUMR(3,NRING) DIMENSION BUFIN(LSAM),CIROLD(LCIRC),CIRNEW(LCIRC) INTEGER :: KLIST(1) DOUBLE PRECISION :: TEMP(MAXRIN,2), & ENER,TOTMIN,TOTMIM,SOLD,SNEW,EAV,EC(NIMA) DIMENSION :: DIST(NIMA),ROT(NIMA) DIMENSION :: ILIST(NIMA),DLIST(NLIST) CHARACTER(LEN=*) :: FINPAT CHARACTER(LEN=MAXNAM) :: FINPIC,FINP INTEGER :: NLET COMMON /MXR/ MAXRI ! DANGER USED IN ANG( FUNCTION al LOGICAL*1 :: CH_ANG,NORM CHARACTER*1 :: MODE,MIRROR C -------------------------------------------- C USED ONLY IN HKMC INTEGER*2 :: IP(NIMA),IQ(NKMAX) DIMENSION :: CIRSEED(LCIRC,NKMAX) DOUBLE PRECISION :: E(NKMAX),ES(NKMAX) C -------------------------------------------- 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), CIRC(LCIRC,NIMA), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN MWANT = NSAM*NROW + LCIRC*NIMA CALL ERRT(46,'AP C, X & CIRC',MWANT) RETURN ENDIF DO K1=1,NIMA CALL FILGET(FINPAT,FINPIC,NLET,ILIST(K1),INTFLAG) IF (INTFLAG .NE. 0) THEN CALL ERRT(18,'AP C ',NE) GOTO 9999 ENDIF MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,INPIC,'O',IFORM, & NSAMT,NROWT,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 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) C NORMALIZE IF REQUESTED IF (NORM) CALL NORMAS(X,-NSAM/2,NSAM/2,-NROW/2,NROW/2, & NUMR,NUMR(1,NRING)) 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 DIST IS USED HERE FOR THE RANDOM CHOOSING OF IMAGES DIST = 0.0 CALL RANDOM_NUMBER(CIID) IMI = MIN(NIMA,MAX(1,INT(CIID*NIMA+0.5))) CIROLD = CIRC(:,IMI) ROT(IMI) = 1.0 DIST(IMI) = 1.0 SOLD = 0.0D0 EAV = EC(IMI) DO KTN=2,NIMA 804 CALL RANDOM_NUMBER(CIID) M = MIN(NIMA, MAX(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 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) IF (MIRROR .EQ. 'M') THEN CALL CROSRMG(CIROLD,CIRC(1,IMI),LCIRC,NRING, & TEMP,TEMP(1,2), & MAXRIN,JACUP,NUMR,TOTMIM,TMT,MODE) IF (TMT .GT. TOT) THEN ROT(IMI) = -TMT SOLD = SOLD+EAV+EC(IMI)-2.0*TOTMIM CALL UPDTM(CIROLD,CIRC(1,IMI),LCIRC,NRING,NUMR,TOT, & MAXRIN,KTN) GOTO 151 ENDIF ENDIF ROT(IMI) = TOT SOLD = SOLD + EAV + EC(IMI) - 2.0 * TOTMIN CALL UPDTC(CIROLD,CIRC(1,IMI),LCIRC,NRING,NUMR,TOT, & MAXRIN,KTN) 151 EAV = ENER(CIROLD,LCIRC,NRING,NUMR,MODE) ENDDO CIRNEW = CIROLD ROT = 0.0 C WRITE(NOUT,*) SOLD*FLOAT(NIMA)/(NIMA-1) C WRITE(NOUT,2001) (ANG(ROT(J),MODE),J=1,NIMA) C ITERATIONS TO GET BETTER APPROXIMATION ITER = 0 901 CONTINUE ITER = ITER+1 CH_ANG = .FALSE. SNEW = 0.0D0 C DO IMI=1,NIMA CALL CROSRNG(CIROLD,CIRC(1,IMI),LCIRC,NRING,TEMP,TEMP(1,2), & MAXRIN,JACUP,NUMR,TOTMIN,TOT,MODE) IF (MIRROR .EQ. 'M') THEN CALL CROSRMG(CIROLD,CIRC(1,IMI),LCIRC,NRING,TEMP, & TEMP(1,2),MAXRIN,JACUP,NUMR,TOTMIM,TMT,MODE) IF (TMT .GT. TOT) THEN IF (ROT(IMI) .NE. -TMT) THEN CH_ANG = .TRUE. ROT(IMI) = -TMT ENDIF CALL UPDTM(CIRNEW,CIRC(1,IMI),LCIRC,NRING,NUMR, & TMT,MAXRIN,IMI) TOTMIN = TOTMIM GOTO 152 ENDIF ENDIF 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) 152 SNEW = SNEW+EAV+EC(IMI)-2.0*TOTMIN DIST(IMI) = EAV+EC(IMI)-2.0*TOTMIN ENDDO WRITE(NOUT,2020) ITER,SNEW 2020 FORMAT(' Iteration #',I3,' Sum of distances=',1PD13.6) IF (SNEW.LE.SOLD .AND. CH_ANG) THEN CIROLD = CIRNEW EAV = ENER(CIROLD,LCIRC,NRING,NUMR,MODE) SOLD = SNEW GOTO 901 ENDIF WRITE(NOUT,2001) (SIGN(ANG(ABS(ROT(J)),MODE),ROT(J)),J=1,NIMA) 2001 FORMAT(8(1X,F8.3)) CALL SEEDS(CIRSEED,CIRC,DIST,NKMAX,LCIRC,IP,NIMA,NOUT) CALL HKMC(CIRSEED,CIRC,CIRNEW,NKMAX,LCIRC,IP,IQ,ES,EC,E, & DIST,ROT,NRING,TEMP,MAXRIN,JACUP,NUMR,MAXIT, & NIMA,MODE,SNEW,NOUT,MIRROR) NMAX = 0 CALL FILSEQP(FINPAT,NLET,KLIST,NMAX,NIXX, & 'OBJECT OUTPUT FILENAME TEMPLATE',IRTFLG) DO II=1,NKMAX CALL FILGET(FINPAT,FINP,NLET,II,IRTFLG) III = 0 NLS = 2 DO IIII=1,NIMA IF (IP(IIII) .EQ. II) THEN III = III + 1 DLIST(1) = III DLIST(2) = ILIST(IIII) IAP = 0 CALL SAVDN1(NDOC,FINP,DLIST,NLS,III-1,IAP) ENDIF ENDDO CLOSE(NDOC) ENDDO I = 0 DO IMI=1,NIMA I = I+1 712 IF (ILIST(I) .EQ. -1) THEN I = I+1 GOTO 712 ENDIF DLIST(1) = IMI DLIST(2) = ILIST(I) DLIST(3) = ANG(ABS(ROT(IMI)),MODE) DLIST(4) = SIGN(DIST(IMI),ROT(IMI)) DLIST(5) = IP(IMI) CALL SAVD(NDOC,DLIST,NLIST,IRTFLG) ENDDO CALL SAVDC CLOSE(NDOC) 9999 IF (ALLOCATED(X)) DEALLOCATE(X) IF (ALLOCATED(CIRC)) DEALLOCATE(CIRC) END C ---------------------- UPDTM ------------------------------ SUBROUTINE UPDTM(CIRC1,CIRC2,LCIRC,NRING,NUMR,TOT,MAXRIN,IS) DIMENSION :: CIRC1(LCIRC),CIRC2(LCIRC) INTEGER :: NUMR(3,NRING),MAXRIN COMPLEX :: C PI2 = 8.0D0*DATAN(1.0D0) c$omp parallel do private(i,j,nsirt,arg,c) DO I=1,NRING NSIRT = NUMR(3,I) CIRC1(NUMR(2,I)) = & (CIRC1(NUMR(2,I))*(IS-1)+CIRC2(NUMR(2,I)))/REAL(IS) CIRC1(NUMR(2,I)+1) = & (CIRC1(NUMR(2,I)+1)*(IS-1)+CIRC2(NUMR(2,I)+1)* & COS(PI2*(TOT-1.0)/2.0 & *REAL(NSIRT)/REAL(MAXRIN)))/REAL(IS) DO J=3,NSIRT,2 ARG = PI2*(TOT-1.0)*REAL(J/2)/REAL(MAXRIN) C = CMPLX(CIRC2(NUMR(2,I)+J-1),-CIRC2(NUMR(2,I)+J))* & CMPLX(COS(ARG),SIN(ARG)) CIRC1(NUMR(2,I)+J-1) = & (CIRC1(NUMR(2,I)+J-1)*(IS-1)+REAL(C))/REAL(IS) CIRC1(NUMR(2,I)+J) = & (CIRC1(NUMR(2,I)+J)*(IS-1)+AIMAG(C))/REAL(IS) ENDDO ENDDO END