C++********************************************************************* C C WRITPRO_N.F SPEEDED UP FEB 2000 ARDEAN LEITH C VERBOSE OUTPUT NOV 2000 ARDEAN LEITH C PUT ANGLES IN HEADER JUN 2001 ARDEAN LEITH C OPFILEC FEB 2003 ARDEAN LEITH C PJ_RRINGS FEB 2005 ARDEAN LEITH C APRINGS_NEW APR 2008 ARDEAN LEITH CC ********************************************************************** 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 WRITPRO_N(PROJ,LUNPROJ,NSAM,NROW,NUMTH,BCKE,NNN, C IPCUBE,NN,RI,ISELECT,NANG,MAXKEY,ANGBUF) C LUNRINGS,MODE,MR,NR,ISKIP,IRTFLG C C PARAMETERS: C NUMTH NUMBER OF OMP THREADS SENT C RI RADIUS SENT C C PURPOSE: COMPUTES A PROJECTION OF A 3D VOLUME ACCORDING TO THE C THREE EULERIAN ANGLES CAN ALSO SAVE REFERENCE RINGS IN C POLAR FORM C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE WRITPRO_N(PROJ,LUNPROJ,NSAM,NROW,NSLICE,NUMTH, & BCKE,NNN,IPCUBE,NN,RI,ISELECT,NANG,MAXKEY,ANGBUF, & LUNRINGS,MODE,MR,NR,ISKIP,IRTFLG) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' INTEGER, DIMENSION(MAXKEY) :: ISELECT REAL, DIMENSION(3,MAXKEY) :: ANGBUF REAL, DIMENSION(NSAM,NROW,NUMTH) :: PROJ REAL, DIMENSION(NNN) :: BCKE REAL, DIMENSION(5,NN) :: IPCUBE REAL, DIMENSION(4) :: BUFOUT CHARACTER(LEN=*) :: MODE CHARACTER(LEN=MAXNAM) :: PROJNAM,PROJPAT,REFNAM LOGICAL :: REFRINGS INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NUMR REAL, ALLOCATABLE, DIMENSION(:) :: WR INTEGER, PARAMETER :: NPLANS = 14 #ifndef SP_32 INTEGER *8 :: FFTW_PLANS(NPLANS) #else INTEGER *4 :: FFTW_PLANS(NPLANS) #endif #ifdef USE_MPI include 'mpif.h' ICOMM = MPI_COMM_WORLD CALL MPI_COMM_RANK(ICOMM, MYPID, IERR) #else MYPID = -1 #endif REFRINGS = (LUNRINGS .GT. 0) IF (REFRINGS) THEN C FIND NUMBER OF REFERENCE-RINGS NRING = ((NR - MR) / ISKIP) + 1 C FILL NUMR WITH RING ADDRESSES ALLOCATE(NUMR(3,NRING), WR(NRING),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(46,'NUMR & WR',4*NRING) GOTO 9999 ENDIF C NUMBR(1,*) IS RING NUMBER NRING = 0 DO I=MR,NR,ISKIP NRING = NRING + 1 NUMR(1,NRING) = I ENDDO C CALCULATES DATA FOR NUMR & LCIRC CALL ALPRBS_Q(NUMR,NRING,LCIRC,MODE) C RINGWE RETURNS WR CALL RINGWE_NEW(WR,NUMR,NRING,NUMR(3,NRING)) IF (MODE .EQ. 'H') WR = WR * 0.5 C INITIALIZE FFTW3 PLANS FOR USE WITHIN OMP || SECTIONS CALL APRINGS_INIT_PLANS(NUMR,NRING, & FFTW_PLANS,NPLANS,IRTFLG) ENDIF C GET PROJPAT TEMPLATE ONLY (NOT ILIST) NMAX = 0 CALL FILSEQP(PROJPAT,NLET,ILIST,NMAX,NIMA, & 'ENTER TEMPLATE FOR 2-D PROJECTION',IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 IF (REFRINGS) THEN C OPEN REF RINGS FILE USING SPIDER FORMAT MAXIM = 0 CALL OPFILEC(0,.TRUE.,REFNAM,LUNRINGS,'U',IFORM, & LCIRC,NANG,1,MAXIM,'REFERENCE RINGS',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 9999 ENDIF IF (VERBOSE .AND. MYPID .LE. 0) THEN WRITE(NOUT,*) ' ' ENDIF IFORM = 1 NSL = 1 IANG = 1 DO WHILE (IANG .LE. NANG) NEEDED = MIN(IANG+NUMTH-1,NANG) C BREAK INTO NUMTH CHUNKS, CALCULATE PROJECTIONS c$omp parallel do private(i,ifile) DO I=IANG,NEEDED IFILE = ISELECT(I) CALL WPRO_N(PROJ(1,1,I-IANG+1),NSAM,NROW,NSLICE,BCKE, & NNN,IPCUBE,NN, & ANGBUF(3,IFILE),ANGBUF(2,IFILE),ANGBUF(1,IFILE),RI) ENDDO C WRITE PROJECTIONS TO OUTPUT FILES DO I=IANG,NEEDED IFILE = ISELECT(I) C CREATE OUTPUT FILENAME CALL FILGET(PROJPAT,PROJNAM,NLET,IFILE,IRTFLG) C OPEN NEXT PROJECTION OUTPUT FILE MAXIM = 0 CALL OPFILEC(0,.FALSE.,PROJNAM,LUNPROJ,'U',IFORM, & NSAM,NROW,NSL,MAXIM,' ',.FALSE.,IRTFLG) CALL WRTVOL(LUNPROJ,NSAM,NROW,1,1, & PROJ(1,1,I-IANG+1),IRTFLG) IF (VERBOSE .AND. MYPID .LE. 0) THEN WRITE(NOUT,90) IFILE, & ANGBUF(1,IFILE),ANGBUF(2,IFILE),ANGBUF(3,IFILE) 90 FORMAT(' PROJECTION:',I6, & ' PSI:',F6.1,' THETA:',F6.1,' PHI:',F6.1) ENDIF C PUT ANGLES IN HEADER BUFOUT(1) = ANGBUF(1,IFILE) BUFOUT(2) = ANGBUF(2,IFILE) BUFOUT(3) = ANGBUF(3,IFILE) BUFOUT(4) = 1.0 CALL LUNSETVALS(LUNPROJ,IAPLOC+1,4,BUFOUT,IRTFLG) CLOSE(LUNPROJ) ENDDO IF (REFRINGS) THEN C CALCULATE & WRITE REFERENCE RINGS CALL PJ_RRINGS(NSAM,NROW, & NRING,LCIRC,NUMR,MODE, & NUMTH,LUNRINGS,WR,FFTW_PLANS, & PROJ,IANG,NEEDED,IRTFLG) ENDIF IANG = IANG + NUMTH ENDDO 9999 IF (REFRINGS) THEN IF (ALLOCATED(NUMR)) DEALLOCATE(NUMR) IF (ALLOCATED(WR)) DEALLOCATE(WR) ENDIF RETURN END C ---------------------- PJ_RRINGS -------------------------- SUBROUTINE PJ_RRINGS(NSAM,NROW, & NRING,LCIRC,NUMR,MODE, & NUMTH,LUNRINGS,WR,FFTW_PLANS, & PROJ,IANG1,IANG2,IRTFLG) CHARACTER(LEN=1) :: MODE C EXTERNAL ARRAYS INTEGER,DIMENSION(3,NRING) :: NUMR REAL,DIMENSION(NSAM,NROW,NUMTH) :: PROJ REAL,DIMENSION(NRING) :: WR C FFTW_PLANS CONTAINS POINTERS TO STRUCTURES INTEGER*8, INTENT(IN) :: FFTW_PLANS(*) C AUTOMATIC ARRAYS REAL,DIMENSION(LCIRC,NUMTH) :: CIRCREF IRTFLG = 0 C CALCULATE CENTER FOR APRINGS CNS2 = NSAM / 2 + 1 CNR2 = NROW / 2 + 1 c$omp parallel do private(IMI,IT) DO IMI=IANG1,IANG2 IT = IMI - IANG1 + 1 C NORMALIZE IMAGE VALUES UNDER THE MASK OVER VARIANCE RANGE C INTERPOLATE TO POLAR COORDINATES, CREATE FOURIER OF: CIRCREF C WEIGHT CIRCREF USING WR CALL APRINGS_ONE_NEW(NSAM,NROW, CNS2,CNR2, & PROJ(1,1,IT),.FALSE., & MODE,NUMR,NRING,LCIRC,WR,FFTW_PLANS, & CIRCREF(1,IT),IRTFLG) ENDDO C SAVE CIRCREF IN FILE OPENED ON LUNRINGS DO IMI=IANG1,IANG2 IT = IMI - IANG1 + 1 CALL WRTLIN(LUNRINGS,CIRCREF(1,IT),LCIRC,IMI) c write(6,*) 'oCIRC(1-100):',CIRCREF(1,1),CIRCREF(100,1) ENDDO END #ifdef NEVER C I TRIED THIS BUT IT WAS A BIT SLOWER !!!!!!!!!!!-------------- C -------------- PJ_RRINGSNEW----------------------------------- SUBROUTINE PJ_RRINGSNEW(NSAM,NROW, & NRING,LCIRC,NUMR,MODE, & NUMTH,LUNRINGS,WR, & PROJ,IANG1,IANG2,IRTFLG) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' CHARACTER(LEN=1) :: MODE C EXTERNAL ARRAYS INTEGER,DIMENSION(3,NRING) :: NUMR REAL,DIMENSION(NSAM,NROW,NUMTH) :: PROJ REAL,DIMENSION(NRING) :: WR C AUTOMATIC ARRAYS REAL,DIMENSION(LCIRC,NUMTH) :: CIRCREF IRTFLG = 0 C PREPARE CIRCULAR RINGS DATA FOR REFERENCE IMAGES c$omp parallel do private(IMI,IT) DO IMI=IANG1,IANG2 IT = IMI - IANG1 + 1 !PROJ & CIRCREF INDEX C NORMALIZE, INTERPOLATE INTO POLAR COORDINATES, FFT, & WEIGHT CALL APRINGS_ONE(NSAM,NROW,WR, & PROJ(1,1,IT),CIRCREF(1,IT), & MODE,NUMR,NRING,LCIRC, IRTFLG) ENDDO C SAVE CIRCREF IN FILE OPENED ON LUNRINGS DO IMI=IANG1,IANG2 IT = IMI - IANG1 + 1 CALL WRTLIN(LUNRINGS,CIRCREF(1,IT),LCIRC,IMI) if (imi.eq.1) THEN write(6,*) 'n(1-100):',CIRCREF(1,1),CIRCREF(100,1),numth endif ENDDO END #endif