C++******************************************************************* C C FFILTS.F USED OPFILE NOV 00 ARDEAN LEITH C ADDED BFACTOR OCT 01 BILL BAXTER C OPFILEC FEB 03 ARDEAN LEITH C GAUSSIAN BUG FEB 04 PP C SAMPLED ADDED MAR 07 C. RENKEN C C ********************************************************************** C=* FROM: SPIDER - MODULAR IMAGE PROCESSING SYSTEM. AUTHOR: J.FRANK * C=* Copyright (C) 1985-2007 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 FFILTS(LUN,NSAM,NROW,NSLICE,NSAMO) C LUN LOGICAL UNIT NUMBER OF FOURIER FILE TO BE FILTERED C NSAM,NROW,NSLICE DIMENSIONS OF FOURIER FILE C NSAMO DIMENSION OF CORRESPONDING REAL-SPACE FILE C C--******************************************************************* SUBROUTINE FFILTS(LUN,LUNO,NSAM,NROW,NSLICE,NSAMO) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' COMMON B(1) EQUIVALENCE (CB,B) COMPLEX CB(1) CHARACTER *(MAXNAM) FILNAM CHARACTER*1 NULL DATA LUNF/27/ NULL = CHAR(0) 951 WRITE(NOUT,1009) 1009 FORMAT & (' 1 - LOW-PASS, 2 - HIGH-PASS',/, & ' 3 - GAUSS LOW-PASS, 4 - GAUSS HIGH-PASS',/, & ' 5 - FERMI LOW-PASS, 6 - FERMI HIGH-PASS',/, & ' 7 - BUTER LOW-PASS, 8 - BUTER HIGH-PASS',/, & ' 9 - REMEZ, 10 - B FACTOR',/, & ' 11 - SAMPLED SPACE') CALL RDPRMI(IOPT,IDUM,NOT_USED,'FILTER TYPE (1-11)') IF (IOPT.LT.1 .OR. IOPT.GT.11) THEN CALL ERRT(102,'NO SUCH FILTER TYPE',IOPT) GOTO 951 ENDIF C B FACTOR FILTER ************************************************* IF (IOPT .EQ. 10) THEN CALL BFACT(LUN,LUNO,NSAM,NROW,NSLICE,NSAMO) RETURN ENDIF C REMEZ FILTER ************************************************* IF (IOPT .EQ. 9) THEN MAXIM = 0 CALL OPFILEC(0,.TRUE.,FILNAM,LUNF,'O',IFORM,NS1,NR1,NSL1, & MAXIM,'FILTER',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (IFORM .GT .0 .OR. & NS1.NE.NSAM .OR. NR1.NE.NROW .OR. NSL1.NE.NSLICE) THEN CALL ERRT(2,'FILT2 ',NE) CLOSE(LUNF) RETURN ENDIF NS2=NSAM/2 DO K=1,NSLICE DO J=1,NROW NR=(K-1)*NROW+J CALL REDLIN(LUN,CB,NSAM,NR ) CALL REDLIN(LUNF,CB(NS2+1),NSAM,NR) DO I=1,NSAM/2 CB(I)=CB(I)*CB(I+NS2) ENDDO CALL WRTLIN(LUNO,CB,NSAM,NR) ENDDO ENDDO CLOSE(LUNF) RETURN ENDIF CC 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) GO TO 5768 ENDIF C FOR SAMPLED SPACE FILTER ******************************* IF (IOPT .EQ. 11)THEN CALL RDPRM(NUMPRJ,NOT_USED,'NUMBER OF PROJECTIONS') GOTO 5768 ENDIF C OTHER FILTERS *********************************************** 952 CALL RDPRM(PARM1,NOT_USED,'FILTER RADIUS') IF (PARM1.LT.0.0.OR.PARM1.GT.0.5) PARM1 = 0.5*PARM1/(NSAMO/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 5768 NS2 = NSAM/2 NR2 = NROW/2 NL2 = NSLICE/2 X1 = FLOAT(NSAMO/2)**2 Y1 = FLOAT(NR2)**2 IF (NSLICE .EQ. 1) THEN Z1 = 1.0 ELSE Z1 = FLOAT(NL2)**2 ENDIF PARM = PARM1**2 DO K=1,NSLICE IZ = K-1 IF (IZ .GT. NL2) IZ = IZ-NSLICE DO J=1,NROW IY = J-1 IF (IY .GT. NR2) IY = IY-NROW NR = J+(K-1)*NROW CALL REDLIN(LUN,B,NSAM,NR ) DO 2 I=1,NS2 IX=I-1 IF (IOPT.EQ.11) GOTO 800 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) CB(I) = CMPLX(0.0,0.0) GOTO 2 C HIGH PASS **************************************************** 200 IF((IX.NE.0 .OR. IY.NE.0 .OR. IZ.NE.0) .AND. & 0.25*(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1) & .LE. PARM) CB(I) = CMPLX(0.0,0.0) 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 CB(I) = CB(I) * EXP(-F) ELSE CB(I) = CMPLX(0.0,0.0) ENDIF GOTO 2 C GAUSSIAN HIGH PASS ******************************************** 400 IF (IX.NE.0 .OR. IY.NE.0 .OR. IZ.NE.0) THEN F=0.125*(FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1)/ & PARM IF (F .LT. 16.0) THEN CB(I) = CB(I) * (1.0-EXP(-F)) ENDIF ENDIF GOTO 2 C FERMI DISTRIBUTION FILTER ************************************* 500 ARG = (0.5*SQRT( & FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1)-PARM1)/TEMP ARG = AMIN1(AMAX1(ARG,-10.0),10.0) IF (IOPT.EQ.6.AND.IX.NE.0.AND.IY.NE.0.AND.IZ.NE.0) GOTO 2 CB(I) = CB(I) * (1.0/(1.0+EXP(ARG))) GOTO 2 C BUTTERWORTH LOWPASS FILTER ********************************** 600 ARG = 0.5*SQRT( & FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1) CB(I) = CB(I) * SQRT(1.0/(1.0+(ARG/PARM1)**ORD)) GOTO 2 C BUTTERWORTH HIGHPASS FILTER ********************************** 700 IF (IX.NE.0 .OR. IY.NE.0 .OR. IZ.NE.0) THEN ARG = 0.5*SQRT( & FLOAT(IX*IX)/X1+FLOAT(IY*IY)/Y1+FLOAT(IZ*IZ)/Z1) CB(I) = CB(I)*(1.0-SQRT(1.0/(1.0+(ARG/PARM1)**ORD))) ENDIF GOTO 2 C SAMPLED SPACE FILTER ***************************************** 800 ARG = SQRT(FLOAT(IX*IX + IY*IY + IZ*IZ)) IF (ARG .EQ. 0.0) ARG = 1.0 F = 3*NUMPRJ*((ARG+0.5)**2 - (ARG-0.5)**2)/ & (4*((ARG+0.5)**3 - (ARG-0.5)**3)) IF (F .GT. 1) F = 1.0 CB(I) = CB(I) * F 2 CONTINUE CALL WRTLIN(LUNO,B,NSAM,NR) ENDDO ENDDO END