C++********************************************************************* C C SAQB_P.F 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 IMAGE_PROCESSING_ROUTINE C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE SAQB_P(BUF,LSD,NSAM,NROW,NIMA,NGRP,JACUP, & CNEW) PARAMETER (NLIST=7) INCLUDE 'CMBLOCK.INC' COMMON /FISPEC/ FINPAT,NLET,FINPIC,DOCFIL,NLETI REAL, ALLOCATABLE, DIMENSION(:,:) :: A,B,C,D REAL, ALLOCATABLE, DIMENSION(:,:,:) ::X DIMENSION CNEW(LSD,NROW) DOUBLE PRECISION SOLD,SNEW,EAV,ENE(NIMA),ENFR_2 DIMENSION DLIST(NLIST),BUF(1024) DIMENSION DIST(NIMA),SHIFT(2,NIMA) LOGICAL CKOT(NIMA) LOGICAL*1 CH_ANG CHARACTER*80 FINPAT,DOCFIL,FINPIC DATA INPIC/77/,NDOC/55/,NDOUT/56/ ALLOCATE (A(LSD,NROW),B(LSD,NROW),C(LSD,NROW), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'AP SA, A,B & C',IER) RETURN ENDIF ALLOCATE (D(LSD,NROW), X(LSD,NROW,NIMA), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'SAQB_P, D & X',IER) DEALLOCATE (A,B,C) RETURN ENDIF NSNR=LSD*NROW SOLD=8.0D0*DATAN(1.0D0)/360.0D0 REWIND NDOC K1=1 KT1=1 KT2=1 IMI=0 778 LERR=-1 CALL UNSAV(DOCFIL,KT1,NDOC,K1,DLIST,4,LERR,KT2) IF (LERR.NE.0) GOTO 788 K1=K1+1 IF (IFIX(DLIST(4)).NE.NGRP) GOTO 778 IMI=IMI+1 C AN(IMI)=DLIST(2) L=DLIST(1) CALL FILGET(FINPAT,FINPIC,NLET,L,INTFLAG) MAXIM = 0 CALL OPFILEC(0,.FALSE.,FINPIC,INPIC,'O',IFORM,NSAMT,NROWT,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) THEN DEALLOCATE (A,B,C,D,X) RETURN ENDIF DO K2=1,NROW CALL REDLIN(INPIC,A(1,K2),NSAM,K2) ENDDO CLOSE(INPIC) CALL RTQ_Q(A,X(1,1,IMI),LSD,NSAM,NROW,DLIST(2)) INS=1 CALL FMRS_2(X(1,1,IMI),NSAM,NROW,INS) IF (INS.EQ.0) THEN CALL ERRT(38,'AP SA',NE) DEALLOCATE (A,B,C,D,X) RETURN ENDIF GOTO 778 788 CONTINUE c$omp parallel do private(imi) DO IMI=1,NIMA ENE(IMI)=ENFR_2(X(1,1,IMI),LSD,NSAM,NROW) ENDDO C BUILD FIRST AVERAGE C TWO ESTIMATION OF INITIAL AVERAGE ARE USED C ONLY ONE !! 11/06/91 DO ICR=1,1 C DIST IS USED HERE FOR THE RANDOM CHOOSING OF IMAGES c$omp parallel do private(imi) DO IMI=1,NIMA DIST(IMI) = 0.0 ENDDO CALL RANDOM_NUMBER(CIID) IMI = MIN0(NIMA,MAX0(1,INT(CIID*NIMA+0.5))) DIST(IMI) = 1.0 CKOT(IMI) = .FALSE. SHIFT(1,IMI) = 0.0 SHIFT(2,IMI) = 0.0 CALL COP(X(1,1,IMI),C,NSNR) 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 ENDDO GOTO 804 810 IMI = I DIST(IMI) = 1.0 LSC = NSAM+2-MOD(NSAM,2) CALL CCRS_2(C,X(1,1,IMI),D, LSC,NSAM,NROW) CALL FINDMX(D,LSD,NSAM,NROW,CMX1,SX1,SY1,JACUP) CALL CR180_2(C,X(1,1,IMI),D,LSD/2,NSAM,NROW) CALL FINDMX(D,LSD,NSAM,NROW,CMX2,SX2,SY2,JACUP) SX2=-SX2 SY2=-SY2 IF(CMX1.GE.CMX2) THEN CKOT(IMI)=.FALSE. SHIFT(1,IMI)=SX1 SHIFT(2,IMI)=SY1 CALL SHFC_2(X(1,1,IMI),A,LSD/2,NSAM,NROW,SX1,SY1) CALL UPDTF(C,A,NSNR,KTN) ELSE CKOT(IMI)=.TRUE. SHIFT(1,IMI)=SX2 SHIFT(2,IMI)=SY2 CALL SH180_2(X(1,1,IMI),A,LSD/2,NSAM,NROW,SX2,SY2) CALL UPDTF(C,A,NSNR,KTN) ENDIF ENDDO SOLD=ENFR_2(C,LSD,NSAM,NROW) WRITE(NOUT,2055) ICR,SOLD 2055 FORMAT(' Initial estimation #',I1, & ' Squared sum=',1PE12.5) C WRITE(NOUT,2001) (CKOT(J),SHIFT(1,J),SHIFT(2,J),J=1,NIMA) 2001 FORMAT(4(1X,L1,2(1X,F8.3))) C IF (ICR.EQ.1) THEN CALL COP(C,CNEW,NSNR) ELSE LSC = NSAM+2-MOD(NSAM,2) CALL CCRS_2(C,CNEW,D, LSC,NSAM,NROW) CALL FINDMX(D,LSD,NSAM,NROW,CMX1,SX1,SY1,JACUP) CALL CR180_2(C,CNEW,D,LSD/2,NSAM,NROW) CALL FINDMX(D,LSD,NSAM,NROW,CMX2,SX2,SY2,JACUP) SX2=-SX2 SY2=-SY2 IF (CMX1.GE.CMX2) THEN CALL SHFM_2(CNEW,LSD/2,NSAM,NROW,SX1,SY1) CALL UPDTF(C,CNEW,NSNR,2) ELSE CALL SH180_2(CNEW,A,LSD/2,NSAM,NROW,SX2,SY2) CALL UPDTF(C,A,NSNR,2) ENDIF SOLD=ENFR_2(C,LSD,NSAM,NROW) WRITE(NOUT,2065) SOLD 2065 FORMAT(' Initial average. Squared sum=',1PE12.5) c$omp parallel do private(imi) DO IMI=1,NIMA SHIFT(1,IMI)=0.0 SHIFT(2,IMI)=0.0 CKOT(IMI)=.TRUE. ENDDO ENDIF ENDDO C ITERATIONS TO GET BETTER APPROXIMATION ITER=0 901 CONTINUE ITER=ITER+1 CH_ANG=.FALSE. DO IMI=1,NIMA C Subtract from the average ... IF (CKOT(IMI)) THEN CALL SH180_2 & (X(1,1,IMI),B,LSD/2,NSAM,NROW,SHIFT(1,IMI),SHIFT(2,IMI)) ELSE CALL SHFC_2 & (X(1,1,IMI),B,LSD/2,NSAM,NROW,SHIFT(1,IMI),SHIFT(2,IMI)) ENDIF CALL SUBAF(C,B,A,NSNR,NIMA) LSC = NSAM+2-MOD(NSAM,2) CALL CCRS_2(A,X(1,1,IMI),D, LSC,NSAM,NROW) CALL FINDMX(D,LSD,NSAM,NROW,CMX1,SX1,SY1,JACUP) CALL CR180_2(A,X(1,1,IMI),D,LSD/2,NSAM,NROW) CALL FINDMX(D,LSD,NSAM,NROW,CMX2,SX2,SY2,JACUP) SX2=-SX2 SY2=-SY2 IF(CMX1.GE.CMX2) THEN IF(CKOT(IMI) & .OR.SHIFT(1,IMI).NE.SX1.OR.SHIFT(2,IMI).NE.SY1) THEN CH_ANG = .TRUE. CKOT(IMI) = .FALSE. SHIFT(1,IMI) = SX1 SHIFT(2,IMI)=SY1 CALL SHFC_2(X(1,1,IMI),B,LSD/2,NSAM,NROW,SX1,SY1) ENDIF CALL UPDTF(CNEW,B,NSNR,IMI) ELSE IF (.NOT. CKOT(IMI) & .OR.SHIFT(1,IMI).NE.SX2.OR.SHIFT(2,IMI).NE.SY2) THEN CH_ANG =.TRUE. CKOT(IMI) =.TRUE. SHIFT(1,IMI) =SX2 SHIFT(2,IMI) =SY2 CALL SH180_2(X(1,1,IMI),B,LSD/2,NSAM,NROW,SX2,SY2) ENDIF CALL UPDTF(CNEW,B,NSNR,IMI) ENDIF ENDDO SNEW = ENFR_2(CNEW,LSD,NSAM,NROW) WRITE(NOUT,2066) ITER,SNEW 2066 FORMAT(' Iteration: ',I4,' Squared sum:',1PE12.5) IF (SNEW .GE. SOLD .AND. CH_ANG) THEN CALL COP(CNEW,C,NSNR) SOLD=SNEW GOTO 901 ENDIF WRITE(NOUT,2001) (CKOT(J),SHIFT(1,J),SHIFT(2,J),J=1,NIMA) REWIND NDOC K1=1 KT1=1 KT2=1 NIM=0 978 LERR=-1 CALL UNSAV(DOCFIL,KT1,NDOC,K1,DLIST,4,LERR,KT2) IF (LERR.NE.0) GOTO 988 K1=K1+1 IF(IFIX(DLIST(4)).NE.NGRP) GOTO 978 NIM=NIM+1 DO I=5,2,-1 DLIST(I)=DLIST(I-1) ENDDO C Calculate distances ... IF (CKOT(NIM)) THEN CALL SH180_2 & (X(1,1,NIM),B,LSD/2,NSAM,NROW,SHIFT(1,NIM),SHIFT(2,NIM)) CALL SUBAF(C,B,A,NSNR,NIMA) CALL CR180_2(A,X(1,1,NIM),D,LSD/2,NSAM,NROW) CALL FINDMX(D,LSD,NSAM,NROW,CMX1,SX2,SY2,JACUP) ELSE CALL SHFC_2 & (X(1,1,NIM),B,LSD/2,NSAM,NROW,SHIFT(1,NIM),SHIFT(2,NIM)) CALL SUBAF(C,B,A,NSNR,NIMA) LSC = NSAM+2-MOD(NSAM,2) CALL CCRS_2(A,X(1,1,NIM),D, LSC,NSAM,NROW) CALL FINDMX(D,LSD,NSAM,NROW,CMX1,SX1,SY1,JACUP) ENDIF EAV=ENFR_2(A,LSD,NSAM,NROW) SNEW=EAV+ENE(NIM)-2*CMX1 DLIST(1)=NIM IF(CKOT(NIM)) DLIST(3)=DLIST(3)+180.0 DLIST(5)=SHIFT(1,NIM) DLIST(6)=SHIFT(2,NIM) DLIST(7)=SNEW CALL SAVD(NDOUT,DLIST,NLIST,IRTFLG) GOTO 978 988 CALL SAVDC CLOSE(NDOUT) DEALLOCATE(A,B,C,D,X) END