C++******************************************************************* C C FQ3_P.F 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 FQ3_P(LUN,LUNF,B,LSD,N2S,N2R,N2L,NSAM,NROW,NSLICE,IOPT) C C PURPOSE: QUICK FILTERING OF THE REAL-SPACE FILE BY IN-CORE FFT C C LUN LOGICAL UNIT NUMBER OF REAL-SPACE FILE TO BE FILTERED C LUNO LOGICAL UNIT NUMBER OF REAL-SPACE OUTPUT FILE C IOPT TYPE OF THE FILTER C B BUFFER C NSAM,NROW,NSLICE DIMENSIONS OF REAL-SPACE FILE C N2S=2*NSAM C N2R=2*NROW C N2L=2*NSLICE C C IMAGE_PROCESSING_ROUTINE C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE FQ3_P & (LUN,LUNO,B,LSD,N2S,N2R,N2L,NSAM,NROW,NSLICE,IOPT) INCLUDE 'CMBLOCK.INC' DIMENSION B(LSD,N2R,N2L) DOUBLE PRECISION AVE C READ IMAGE DO K=1,NSLICE DO I=1,NROW NR=(K-1)*NROW+I CALL REDLIN(LUN,B(1,I,K),NSAM,NR) ENDDO ENDDO C BORDER PADDING IF (N2S.NE.NSAM.AND.N2R.NE.NROW.AND.N2L.NE.NSLICE) THEN AVE=(SUM(B(1:NSAM,1:NROW,1))+SUM(B(1:NSAM,1:NROW,NSLICE)) & +SUM(B(1:NSAM,1,2:NSLICE-1))+SUM(B(1:NSAM,NROW,2:NSLICE-1)) & +SUM(B(1,2:NROW-1,2:NSLICE-1))+SUM(B(NSAM,2:NROW-1,2:NSLICE-1))) & /REAL(4*(NSAM+NROW+NSLICE)-16) c$omp parallel do private(i,j,k),reduction(+:ave) DO K=1,NSLICE DO J=1,NROW DO I=1,NSAM AVE=AVE+B(I,J,K) ENDDO ENDDO ENDDO AVE = AVE/FLOAT(NSAM)/FLOAT(NROW)/FLOAT(NSLICE) c$omp parallel do private(i,j,k) DO K=1,NSLICE DO J=1,N2R DO I=NSAM+1,N2S B(I,J,K)=AVE ENDDO ENDDO DO J=NROW+1,N2R DO I=1,NSAM B(I,J,K)=AVE ENDDO ENDDO ENDDO c$omp parallel do private(i,j,k) DO K=NSLICE+1,N2L DO J=1,N2R DO I=1,N2S B(I,J,K)=AVE ENDDO ENDDO ENDDO ENDIF INV = 1 CALL FMRS_3(B,N2S,N2R,N2L,INV) IF (INV .EQ. 0)THEN IOPT = -1 RETURN ENDIF C BUTTERWORTH FILTER*********************** IF (IOPT.EQ.7 .OR. IOPT.EQ.8) THEN EPS = 0.882 AA = 10.624 CALL RDPRM2(FP,FS,NOT_USED, & 'PASS-BAND FREQUENCY & STOP-BAND FREQUENCY') ORD = 2.*ALOG10(EPS/SQRT(AA**2-1.0)) ORD = ORD/ALOG10(FP/FS) PARM1 = FP/(EPS)**(2./ORD) ELSE CALL RDPRM(PARM1,NOT_USED,'FILTER RADIUS') IF (PARM1.LT.0.0.OR.PARM1.GT.0.5) PARM1=0.5*PARM1/(NSAM/2) C FERMI DISTRIBUTION FILTER ******************** IF (IOPT.EQ.5 .OR. IOPT.EQ.6) THEN CALL RDPRM(TEMP,NOT_USED,'TEMPERATURE(0=CUTOFF)') C EXPONENTIAL FOR HIGH-PASS OPTION IF(IOPT.EQ.6) TEMP=-TEMP ENDIF ENDIF NR2 = N2R/2 NL2 = N2L/2 X1 = FLOAT(N2S/2)**2 Y1 = FLOAT(NR2)**2 Z1 = FLOAT(NL2)**2 PARM = PARM1**2 C KEEP ZERO TERM FOR HIGH PASS OPTIONS AVG = B(1,1,1) c$omp parallel do private(i,j,k,ix,iy,iz,f) DO K=1,N2L IZ = K-1 IF (IZ.GT.NL2) IZ=IZ-N2L DO J=1,N2R IY = J-1 IF (IY.GT.NR2) IY=IY-N2R DO 2 I=1,LSD,2 IX = (I-1)/2 GOTO(100,200,300,400,500,500,600,700),IOPT C LOWPASS ******************************** 100 IF (0.25*(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1) & .GT. PARM) THEN B(I,J,K) = 0.0 B(I+1,J,K) = 0.0 ENDIF GOTO 2 C HIGH PASS ****************** 200 IF( 0.25*(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1) & .LE. PARM) THEN B(I,J,K) = 0.0 B(I+1,J,K) = 0.0 ENDIF GOTO 2 C GAUSSIAN LOW PASS *************************** 300 F=0.125*(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1)/PARM IF (F.LT.16.0) THEN F = EXP(-F) B(I,J,K) = B(I,J,K)*F B(I+1,J,K) = B(I+1,J,K)*F ELSE B(I,J,K) = 0.0 B(I+1,J,K) = 0.0 ENDIF GOTO 2 C GAUSSIAN HIGH PASS ****************************** 400 F=0.125* (FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1)/PARM IF (F.LT.16.0) THEN F = (1.0-EXP(-F)) B(I,J,K) = B(I,J,K)*F B(I+1,J,K) = B(I+1,J,K)*F ENDIF GOTO 2 C FERMI DISTRIBUTION FILTER ******************** 500 F=(0.5*SQRT( & FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1)-PARM1)/TEMP F = AMIN1(AMAX1(F,-10.0),10.0) F = (1.0/(1.0+EXP(F))) B(I,J,K) = B(I,J,K)*F B(I+1,J,K) = B(I+1,J,K)*F GO TO 2 C BUTTERWORTH LOWPASS FILTER ********************** 600 F=0.5*SQRT( & FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1) F = SQRT(1.0/(1.0+(F/PARM1)**ORD)) B(I,J,K) = B(I,J,K)*F B(I+1,J,K) = B(I+1,J,K)*F GO TO 2 C BUTTERWORTH HIGIPASS FILTER ********************* 700 F = 0.5*SQRT( & FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1) F = (1.0-SQRT(1.0/(1.0+(F/PARM1)**ORD))) B(I,J,K) = B(I,J,K)*F B(I+1,J,K) = B(I+1,J,K)*F 2 CONTINUE ENDDO ENDDO C RESTORE ZERO TERM FOR HIGH PASS OPTIONS IF (IOPT.EQ.2 .OR. IOPT.EQ.4 .OR. IOPT.EQ.6 .OR. IOPT.EQ.8) & B(1,1,1) = AVG C WRITE IMAGE INV = -1 CALL FMRS_3(B,N2S,N2R,N2L,INV) IF (INV.EQ.0)THEN CALL ERRT(38,'FQ',NE) RETURN ENDIF DO K=1,NSLICE DO I=1,NROW NR = (K-1)*NROW+I CALL WRTLIN(LUNO,B(1,I,K),NSAM,NR) ENDDO ENDDO END