C++********************************************************************* C C GALI_P.F DYNAMIC MEMORY ALLOCATION SEPT 2000 BIMAL RATH C REPLACED 'INCORE' SEPT 2001 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 GALI_P C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE GALI_P(LSD,NSAM,NROW,NSI,ILIST,NIMA,MODE,JACUP, & LCIRC,NUMR,NRING,MAXRIN,BLOB,USEBLOB,NR,NOUT) REAL, ALLOCATABLE, DIMENSION(:,:) :: A, B, C, D, E REAL, ALLOCATABLE, DIMENSION(:,:) :: REFER, REFERN, REFERTMP REAL, ALLOCATABLE, DIMENSION(:,:) :: PARA REAL, ALLOCATABLE, DIMENSION(:) :: A_CIRC, REFER_CIRC REAL, ALLOCATABLE, DIMENSION(:) :: QIMAGES DIMENSION ILIST(NIMA), BLOB(LSD,NROW) LOGICAL CHANGE, NOCHANGE, USEBLOB, IMAGE(NIMA), INCORE INTEGER NUMR(3,NRING) C TEMP IS AN AUTOMATIC ARRAY DOUBLE PRECISION TEMP(MAXRIN+2,2),FNRM,DNRM CHARACTER*1 MODE DATA ZERO/0.0/ DATA INPIC/77/ NSNR = LSD * NROW ALLOCATE(A(LSD,NROW), B(LSD,NROW),C(LSD,NROW), D(LSD,NROW), & E(LSD,NROW), REFER(LSD,NROW), REFERN(LSD,NROW), & REFERTMP(LSD,NROW), PARA(3,NIMA), A_CIRC(LCIRC), & REFER_CIRC(LCIRC), STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'AP SR, A,B,C,D,....',IER) GOTO 9999 ENDIF INCORE = .FALSE. ISIZE = NSAM * NROW * NIMA ALLOCATE(QIMAGES(ISIZE), STAT=IRTFLG) IF (IRTFLG .EQ. 0) THEN WRITE(NOUT,*) ' -- IMAGES IN_CORE VERSION -- ' C LOAD IMAGE DATA INTO QIMAGES TO SPEED THINGS UP DO IMI = 1,NIMA CALL GETIMA_Q(QIMAGES(NSAM*NROW*(IMI-1)+1),NSAM,NSAM,NROW, & IMI,ILIST,NIMA,QIMAGES,INPIC,INCORE) ENDDO INCORE = .TRUE. ELSE WRITE(NOUT,*) ' -- IMAGES FROM DISK VERSION -- ' ENDIF C PREPARE 'BLOB' TO CENTER INPUT IMAGES .... IF (USEBLOB) CALL BLOB_Q(BLOB,LSD,NSAM,NROW,NR) INS = 1 CALL FMRS_2(BLOB,NSAM,NROW,INS) IF (INS .EQ. 0) THEN CALL ERRT(38,'AP SR',NE) GOTO 9999 ENDIF C BUILD FIRST AVERAGE C IMAGE IS USED HERE FOR THE RANDOM SELECTION OF IMAGES C INITIALIZE ALL IMAGE() TO .FALSE. IMAGE = .FALSE. CALL RANDOM_NUMBER(CIID) IMI = MIN0(NIMA,MAX0(1,INT(CIID*NIMA+0.5))) IMAGE(IMI) = .TRUE. PARA(1,IMI) = 0.0 PARA(2,IMI) = 0.0 PARA(3,IMI) = 0.0 C GET IMAGE DATA CALL GETIMA_Q(A,LSD,NSAM,NROW,IMI,ILIST,NIMA, & QIMAGES,INPIC,INCORE) C COPIES IMAGE DATA FROM A INTO B CALL COP(A,B,NSNR) C FT ON IMAGE DATA IN ARRAY B INS = 1 CALL FMRS_2(B,NSAM,NROW,INS) C BLOB - BLOB; B-CURRENT IMAGE; C-CCF. LSC = NSAM+2-MOD(NSAM,2) CALL CCRS_2(BLOB,B,C, LSC,NSAM,NROW) CALL FINDMX_Q(C,LSD,NSAM,NROW,NSI,CMX1,SX1,SY1) ISX1 = SX1 ISY1 = SY1 PARA(2,IMI) = ISX1 PARA(3,IMI) = ISY1 CALL SHFI_Q(A,REFERN,LSD,NSAM,NROW,ISX1,ISY1) C PRINT *,' FIRST IMAGE WITH BLOB',IMI,ISX1,ISY1 CALL COP(REFERN,REFER,NSNR) CALL ALRQ_Q(REFER,LSD,NSAM,NROW,NUMR,REFER_CIRC,LCIRC, & NRING,MODE,IPIC) CALL FOURING_Q(REFER_CIRC,LCIRC,NUMR,NRING,TEMP,MODE) INS = 1 CALL FMRS_2(REFER,NSAM,NROW,INS) C GO THROUGH ALL THE REMAINING IMAGES. CHANGE = .FALSE. DO KTN = 2,NIMA DO CALL RANDOM_NUMBER(CIID) M = MIN0(NIMA,MAX0(1,INT(CIID*(NIMA-KTN+1)+0.5))) IMI = 0 DO I=1,NIMA IF (.NOT.IMAGE(I)) THEN IMI = IMI + 1 IF (IMI .EQ. M) THEN IMI = I GOTO 801 ENDIF ENDIF ENDDO ENDDO 801 IMAGE(IMI) = .TRUE. C GET IMAGE DATA CALL GETIMA_Q(A,LSD,NSAM,NROW,IMI,ILIST,NIMA, & QIMAGES,INPIC,INCORE) C COPY IMAGE DATA ARRAY A INTO ARRAY B CALL COP(A,B,NSNR) INS = 1 C FT ON IMAGE DATA IN ARRAY B CALL FMRS_2(B,NSAM,NROW,INS) C BLOB - BLOB; B-CURRENT IMAGE; C-CCF. LSC = NSAM+2-MOD(NSAM,2) CALL CCRS_2(BLOB,B,C, LSC,NSAM,NROW) CALL FINDMX_Q(C,LSD,NSAM,NROW,NSI,CMX1,SX1,SY1) ISX1 = SX1 ISY1 = SY1 PARA(1,IMI) = 0.0 PARA(2,IMI) = ISX1 PARA(3,IMI) = ISY1 CALL SHFI_Q(A,B,LSD,NSAM,NROW,ISX1,ISY1) C PRINT *,' NEXT IMAGE WITH BLOB',IMI,ISX1,ISY1 C RUN ROTATION-SHIFT ALIGNMENT C INPUT: THREE REAL IMAGES, A, B AND REFERN C OUTPUT: B - REAL ALIGNED IMAGE, PARA - UDATED PARAMETERS C SCRATCH: D C NOCHANGE=.TRUE. - NO CHANGE IN THE POSITION OF THE IMAGE CALL ALROSI_Q(A,B,D,C,REFER,LSD,NSAM,NROW,NSI, & PARA(1,IMI),NOCHANGE, & A_CIRC,REFER_CIRC,LCIRC,JACUP,NUMR,NRING,MAXRIN, & TEMP,MODE,KTN) IF (.NOT. NOCHANGE) CHANGE = .TRUE. C ADD IMAGE TO THE REFERENCE C CALL UPDTF(REFERN,B,NSNR,KTN) CALL UPDTF_R(REFERN,B,NSAM,NROW,LSD,KTN) ENDDO C CORRECT CENTER IF (USEBLOB) THEN CALL CENT_Q(REFERN,LSD,NSAM,NROW,XS,YS) XS = -(XS-NSAM/2-1) YS = -(YS-NROW/2-1) ELSE C BLOB - BLOB; REFERN-CURRENT AVERAGE; D-CCF. CALL COP(REFERN,REFER,NSNR) INS = +1 CALL FMRS_2(REFER,NSAM,NROW,INS) LSC = NSAM+2-MOD(NSAM,2) CALL CCRS_2(BLOB,REFER,D, LSC,NSAM,NROW) CALL FINDMX_Q(D,LSD,NSAM,NROW,NSI,CMX1,XS,YS) ENDIF C PRINT *, ' CENTER CORRECTION ',XS,YS CALL RTQS_Q(REFERN,REFER,LSD,NSAM,NROW,ZERO,XS,YS) DO IMI=1,NIMA PARA(2,IMI) = PARA(2,IMI)+XS PARA(3,IMI) = PARA(3,IMI)+YS ENDDO CALL OUTIM_Q(REFER,LSD,NSAM,NROW,1) CALL OUTPR(PARA,NIMA,1) TEMP(1,1) = FNRM(REFER,NSNR) I = 1 WRITE(NOUT,5091) I,TEMP(1,1) 5091 FORMAT(' ITERATION #',I5,' SUM OF SQUARES=',1PD12.5) C IF NO CHANGES TERMINATE IF (.NOT. CHANGE) GOTO 9999 C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C REFINE THE ALIGNMENT C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ITER = 0 DNRM = 0.0D0 C ITERATE THE REFINEMENT DO CHANGE = .FALSE. ITER = ITER+1 c$omp parallel do private(j,i) DO J=1,NROW DO I=1,LSD REFERN(I,J) = 0.0 ENDDO ENDDO c$omp end parallel do c$omp parallel do private(imi) DO IMI=1,NIMA IMAGE(IMI) = .FALSE. ENDDO c$omp end parallel do DO KTN=1,NIMA DO CALL RANDOM_NUMBER(CIID) M = MIN0(NIMA,MAX0(1,INT(CIID*(NIMA-KTN+1)+0.5))) IMI = 0 DO I=1,NIMA IF (.NOT. IMAGE(I)) THEN IMI = IMI+1 IF (IMI.EQ.M) THEN IMI = I GOTO 901 ENDIF ENDIF ENDDO ENDDO 901 IMAGE(IMI) = .TRUE. C PRINT *,' IMAGE # ',IMI C GET IMAGE DATA IN ARRAY A CALL GETIMA_Q(A,LSD,NSAM,NROW,IMI,ILIST,NIMA, & QIMAGES,INPIC,INCORE) CALL RTQS_Q(A,C,LSD,NSAM,NROW,PARA(1,IMI), & PARA(2,IMI),PARA(3,IMI)) C SUBTRACT THE IMAGE FROM THE REFERENCE CALL SUBAF(REFER,C,REFERTMP,NSNR,NIMA) CALL COP(REFERTMP,E,NSNR) C RUN ROTATION-SHIFT ALIGNMENT C INPUT: TWO REAL IMAGES, B AND REFER C OUTPUT: C - REAL ALIGNED IMAGE, PARA - UDATED PARAMETERS C NOCHANGE=.TRUE. - NO CHANGE IN THE POSITION OF THE IMAGE CALL ALROSF_Q(A,C,D,B,REFERTMP,LSD,NSAM,NROW,NSI, & PARA(1,IMI),NOCHANGE, & A_CIRC,REFER_CIRC,LCIRC,JACUP,NUMR,NRING,MAXRIN,TEMP,MODE) IF (.NOT. NOCHANGE) THEN CHANGE = .TRUE. C ADD IMAGE TO THE REFERENCE CALL UPDF(E,C,REFER,NSNR,NIMA) ENDIF C CALL UPDTF(REFERN,C,NSNR,KTN) CALL UPDTF_R(REFERN,C,NSAM,NROW,LSD,KTN) ENDDO C CORRECT CENTER IF (USEBLOB) THEN CALL CENT_Q(REFERN,LSD,NSAM,NROW,XS,YS) XS = -(XS-NSAM/2-1) YS = -(YS-NROW/2-1) ELSE C BLOB - BLOB; REFERN-CURRENT AVERAGE; D-CCF. CALL COP(REFERN,REFER,NSNR) INS = +1 CALL FMRS_2(REFER,NSAM,NROW,INS) LSC = NSAM+2-MOD(NSAM,2) CALL CCRS_2(BLOB,REFER,D, LSC,NSAM,NROW) CALL FINDMX_Q(D,LSD,NSAM,NROW,NSI,CMX1,XS,YS) ENDIF CALL RTQS_Q(REFERN,REFER,LSD,NSAM,NROW,ZERO,XS,YS) c$omp parallel do private(imi) DO IMI=1,NIMA PARA(2,IMI) = PARA(2,IMI) + XS PARA(3,IMI) = PARA(3,IMI) + YS ENDDO c$omp end parallel do CALL OUTIM_Q(REFER,LSD,NSAM,NROW,ITER+1) CALL OUTPR(PARA,NIMA,ITER+1) TEMP(1,1) = FNRM(REFER,NSNR) WRITE(NOUT,5091) ITER+1,TEMP(1,1) IF (TEMP(1,1) .GE. DNRM) THEN DNRM = TEMP(1,1) IF (.NOT.CHANGE) GOTO 9999 ELSE GOTO 9999 ENDIF C END OF INFINITE ITERATIONS END DO 9999 IF (ALLOCATED(A)) DEALLOCATE(A) IF (ALLOCATED(B)) DEALLOCATE(B) IF (ALLOCATED(C)) DEALLOCATE(C) IF (ALLOCATED(D)) DEALLOCATE(D) IF (ALLOCATED(E)) DEALLOCATE(E) IF (ALLOCATED(REFER)) DEALLOCATE(REFER) IF (ALLOCATED(REFERN)) DEALLOCATE(REFERN) IF (ALLOCATED(REFERTMP)) DEALLOCATE(REFERTMP) IF (ALLOCATED(PARA)) DEALLOCATE(PARA) IF (ALLOCATED(A_CIRC)) DEALLOCATE(A_CIRC) IF (ALLOCATED(REFER_CIRC)) DEALLOCATE(REFER_CIRC) IF (ALLOCATED(QIMAGES)) DEALLOCATE(QIMAGES) END