C++********************************************************************* C C SAQB_F.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_F(BUF,LSD,NSAM,NROW,NIMA,NGRP,JACUP,CNEW) PARAMETER (NLIST=7) INCLUDE 'CMBLOCK.INC' CHARACTER*80 FINPAT,DOCFIL,FINPIC 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_POS DATA IOUTMP/77/ DATA INPIC/77/,NDOC/55/,NDOUT/56/ ALLOCATE (A(LSD,NROW), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'AP SA, A',IER) RETURN ENDIF ALLOCATE (B(LSD,NROW), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'SAQB_P, B',IER) DEALLOCATE (A) RETURN ENDIF ALLOCATE (C(LSD,NROW), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'SAQB_P, C',IER) DEALLOCATE (A,B) RETURN ENDIF ALLOCATE (D(LSD,NROW), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'SAQB_P, D',IER) DEALLOCATE (A,B,C) RETURN ENDIF ALLOCATE (X(LSD,NROW,NIMA), STAT=IRTFLG) IF (IRTFLG.NE.0) THEN CALL ERRT(46,'SAQB_P, X',IER) DEALLOCATE (A,B,C,D) 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',ITYPE,NSAMT,NROWT,NSL, & MAXIM,'DUMMY',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(4,'SAQB_F',NE) 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) 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) 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) (SHIFT(1,J),SHIFT(2,J),J=1,NIMA) 2001 FORMAT(4(2X,2(1X,F8.3))) 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 SHFM_2(CNEW,LSD/2,NSAM,NROW,SX1,SY1) CALL UPDTF(C,CNEW,NSNR,2) 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 ENDDO ENDIF ENDDO C ITERATIONS TO GET BETTER APPROXIMATION ITER=0 901 CONTINUE ITER=ITER+1 CH_POS=.FALSE. DO IMI=1,NIMA C SUBTRACT FROM THE AVERAGE ... CALL SHFC_2 & (X(1,1,IMI),B,LSD/2,NSAM,NROW,SHIFT(1,IMI),SHIFT(2,IMI)) 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) IF(SHIFT(1,IMI).NE.SX1.OR.SHIFT(2,IMI).NE.SY1) THEN CH_POS=.TRUE. 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) 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_POS) THEN CALL COP(CNEW,C,NSNR) SOLD=SNEW GOTO 901 ENDIF WRITE(NOUT,2001) (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 ... 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) EAV=ENFR_2(A,LSD,NSAM,NROW) SNEW=EAV+ENE(NIM)-2*CMX1 DLIST(1)=NIM 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