C++******************************************************************** C C MASK.F ADDED INNER BI SPETP 99 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 MASK(LUNI,LUNO,NSAM,NROW,B) C C PARAMETERS: C LUNI LOGICAL UNIT NUMBER OF INPUT FILE (SENT) C LUNO LOGICAL UNIT NUMBER OF OUTPUT FILE (SENT) C NSAM,NROW,NSLICE DIMENSIONS OF FILE (SENT) C B AVERAGE OF INPUT FILE (SENT / ALTREED) C C--******************************************************************** SUBROUTINE MASK(LUNI,LUNO,NSAM,NROW,NSLICE,B) INCLUDE 'CMBLOCK.INC' DIMENSION BUF(NSAM) PARAMETER (QUADPI = 3.141592653589793238462643383279502884197) CHARACTER * 1 MODE,NULL,ANS DOUBLE PRECISION DAV,AVC,AVCI IF (FCHAR(4:4) .EQ. 'L') THEN C MASK A LINE ACROSS IMAGE IF (NSLICE .GT. 1) THEN C 3-D FILE CALL RDPRMI(NSLICER,NDUM,NOT_USED,'ENTER SLICE NUMBER') ELSE NSLICER = 1 ENDIF CALL RDPRMI(NROWR,NDUM,NOT_USED, & 'ENTER LINE NUMBER TO BE MASKED') CALL RDPRM(B,NOT_USED,'ENTER BACKGROUND') IRECT = (NSLICER - 1) * NROW + NROWR DO I = 1,NROW*NSLICE CALL REDLIN(LUNI,BUF,NSAM,I) IF (I .EQ. IRECT) THEN C MASK THIS LINE BUF = B ENDIF CALL WRTLIN(LUNO,BUF,NSAM,I) ENDDO RETURN ELSE RAD = 0.0 RADI = 0.0 CALL RDPRM2S(RAD,RADI,NOT_USED,'OUTER & INNER RADII',IRTFLG) IF (IRTFLG .NE. 0) RETURN IF (RAD .EQ. 0.0) THEN C FOR INNER ONLY MASKING, SET OUTER RADIUS TO HUGE VALUE RAD = NSAM * 20 ENDIF IF (RAD.LT.0.0 .OR. RADI.LT.0.0 .OR. RAD.LT.RADI) THEN CALL ERRT(14,'MA',IER) RETURN ENDIF ENDIF CALL RDPRMC(MODE,NCHAR,.TRUE., & '(D)isk, (C)osine, (G)aussian edge, or (T)rue Gaussian', & NULL,IRTFLG) IF (IRTFLG .NE. 0) RETURN C BACKGROUND CHOICES: C DAV BACKGROUND IS SET EQUAL TO THE AVERAGE OF THE C IMAGE BEFORE MASKING C PREC AV BACKGROUND IS SET EQUAL TO THE AVERAGE OF THE C IMAGE AREA PASSED BY THE MASK C CIRCUMF BACKGROUND IS SET EQUAL TO THE AVERAGE OF THE C PIXELS ALONG THE MASK'S CIRCUMFERENCE C EXTERNAL BACKGROUND IS SET TO A VALUE SUPPLIED EXTERNALLY CALL RDPRMC(ANS,NCHAR,.TRUE., & '(A)V, (P)REC AV, (C)IRCUMF, OR (E)XTERNAL',NULL,IRTFLG) IF (ANS .EQ. 'E') CALL RDPRM(B,NOT_USED,'ENTER BACKGROUND') CALL RDPRMI(NSAMR,NROWR,NOT_USED, & 'MASK CENTER COORDINATES (NSAM, NROW)') IF (NSLICE .GT. 1) THEN CALL RDPRMI(NSLICER,NDUMP,NOT_USED, & 'MASK CENTER COORDINATE (NSLICE)') ELSE NSLICER = 1 ENDIF C ALTHOUGH THE MASK CENTER COORDINATES MAY BE OUTSIDE THE IMAGE C ITSELF, CHECK TO BE SURE THAT THE RADIUS IS LARGE ENOUGH SO C THAT THE MASK WILL AT LEAST HAVE AN EFFECT ON THE IMAGE. C (APPLIES ONLY TO EXCLUSIVE OUTSIDE MASKING) IF (RAD .GT. 0) THEN C POSITIVE OUTER RADIUS GIVEN IF (((NSAMR + RAD) .LT. 0.0) .OR. & ((NSAMR - RAD - NSAM) .GE. 0.0) .OR. & ((NROWR + RAD) .LT. 0.0) .OR. & ((NROWR - RAD - NROW) .GE. 0.0) .OR. & ((NSLICER + RAD) .LT. 0.0) .OR. & ((NSLICER - RAD-NSLICE) .GE. 0.0)) THEN WRITE(NOUT,61) 61 FORMAT(' *** MASK OUTSIDE IMAGE; NO ACTION TAKEN') RETURN ENDIF ELSEIF (RAD .LE. 0.0) THEN C INNER ONLY MASKING, SET OUTER RADIUS TO HUGE VALUE RAD = NSAM * 20 ENDIF RAD2 = RAD**2 RADI2 = RADI**2 RAD2P = (RAD+1)**2 RADI2P = (RADI-1)**2 C RAD2P IS USED TO DEFINE A NON-EMPTY SET OF POINTS C ALONG THE MASK'S CIRCUMFERENCE IF (ANS .EQ. 'P' .OR. ANS .EQ. 'C') THEN C "P" OR "C" OPTION DAV = 0.0 NAV = 0.0 AVC = 0.0 NAVC = 0.0 AVCI = 0.0 NAVCI = 0.0 DO J = 1,NSLICE FI2 = FLOAT(J-NSLICER)**2 DO I = 1,NROW FI1 = FLOAT(I-NROWR)**2+FI2 CALL REDLIN(LUNI,BUF,NSAM,I+(J-1)*NROW) DO K = 1,NSAM CRAD2 = FI1 + FLOAT(K-NSAMR)**2 IF (CRAD2 .LE. RAD2 .AND. CRAD2 .GE. RADI2) THEN DAV = DAV + BUF(K) NAV = NAV + 1 ENDIF C PROVISION FOR CIRCUMFERENCE OPTION IF (ANS .EQ. 'C') THEN IF (CRAD2 .GE. RAD2 .AND. CRAD2 .LE. RAD2P) THEN C FIND OUTER CIRCUMFERENCE AVC = AVC + BUF(K) NAVC = NAVC + 1 ENDIF IF (RADI2 .GT. 0.0 .AND. & CRAD2 .GE. RADI2P .AND. CRAD2 .LE. RADI2) THEN C FIND INNER CIRCUMFERENCE AVCI = AVCI + BUF(K) NAVCI = NAVCI + 1 ENDIF ENDIF ENDDO ENDDO ENDDO DAV = DAV / NAV B = DAV IF (ANS .EQ. 'C') THEN IF (NAVC .GT. 0) THEN B = AVC / NAVC WRITE(NOUT,21) B 21 FORMAT(' AVERAGE ALONG OUTER CIRCUMFERENCE: ',G12.4) ENDIF IF (NAVCI .GT. 0) THEN BI = AVCI / NAVCI WRITE(NOUT,22) BI 22 FORMAT(' AVERAGE ALONG INNER CIRCUMFERENCE: ',G12.4) ENDIF ENDIF ENDIF C COSINE, GAUSSIAN OR STRAIGHT CUTOFF (DISK) ... C MASKS IN ONE DIRECTION IF (FCHAR(4:4) .EQ. 'X') THEN SWITCHZ = 0.0 SWITCHY = 0.0 SWITCHX = 1.0 ELSEIF (FCHAR(4:4) .EQ. 'Y') THEN SWITCHZ = 0.0 SWITCHY = 1.0 SWITCHX = 0.0 ELSEIF (FCHAR(4:4) .EQ. 'Z') THEN SWITCHZ = 1.0 SWITCHY = 0.0 SWITCHX = 0.0 ELSE SWITCHZ = 1.0 SWITCHY = 1.0 SWITCHX = 1.0 ENDIF IF (MODE .EQ. 'C' ) THEN C COSINE EDGE MASKING CALL RDPRM(HW,NOT_USED,'WIDTH') DO J = 1,NSLICE FI2 = FLOAT(J-NSLICER)**2*SWITCHZ DO I = 1,NROW FI1 = FLOAT(I-NROWR)**2 *SWITCHY+FI2 CALL REDLIN(LUNI,BUF,NSAM,I+(J-1)*NROW) DO K = 1,NSAM CRAD2 = FI1 + SWITCHX*FLOAT(K-NSAMR)**2 IF (CRAD2 .GT. RAD2) THEN SRAD2 = SQRT(CRAD2) WGH = & (1.0+COS(QUADPI*AMIN1(1.0,ABS(SRAD2-RAD)/HW)))*0.5 BUF(K) = B+(BUF(K)-B)*WGH ELSEIF (CRAD2 .LT. RADI2) THEN SRAD2 = SQRT(CRAD2) WGH = & (1.0+COS(QUADPI*AMIN1(1.0,ABS(SRAD2-RADI)/HW)))*0.5 BUF(K) = B+(BUF(K)-B)*WGH ENDIF ENDDO CALL WRTLIN(LUNO,BUF,NSAM,I+(J-1)*NROW) ENDDO ENDDO ELSE IF (MODE .EQ. 'G' ) THEN C GAUSSIAN EDGE MASKING CALL RDPRM(HW,NOT_USED,'HALFWIDTH') HW = -1./(HW**2) DO J = 1,NSLICE FI2 = FLOAT(J-NSLICER)**2*SWITCHZ DO I = 1,NROW FI1 = FLOAT(I-NROWR)**2 *SWITCHY + FI2 CALL REDLIN(LUNI,BUF,NSAM,I+(J-1)*NROW) DO K = 1,NSAM CRAD2 = FI1 + SWITCHX*FLOAT(K-NSAMR)**2 IF (CRAD2.GT.RAD2) THEN SRAD2 = SQRT(CRAD2) WGH = HW*(SRAD2-RAD)**2 IF (WGH .LT. -50.0) THEN WGH = 0.0 ELSE WGH = EXP(WGH) ENDIF BUF(K) = B+(BUF(K)-B)*WGH ELSEIF (CRAD2 .LT. RADI2) THEN SRAD2 = SQRT(CRAD2) WGH=HW*(SRAD2-RADI)**2 IF (WGH.LT.-50.0) THEN WGH = 0.0 ELSE WGH = EXP(WGH) ENDIF BUF(K) = B+(BUF(K)-B)*WGH ENDIF ENDDO CALL WRTLIN(LUNO,BUF,NSAM,I+(J-1)*NROW) ENDDO ENDDO ELSE IF (MODE .EQ. 'T' ) THEN C TRUE GAUSSIAN MASKING IF (RADI .NE. 0.0) THEN WRITE(NOUT,*) ' NO INNER MASKING FOR TRUE GAUSSIAN MASK' CALL ERRT(14,'MA',IER) RETURN ENDIF HW = -1.0 / (RAD**2) DO J = 1,NSLICE FI2 = SWITCHZ*FLOAT(J-NSLICER)**2 DO I = 1,NROW FI1 = SWITCHY*FLOAT(I-NROWR)**2+FI2 CALL REDLIN(LUNI,BUF,NSAM,I+(J-1)*NROW) DO K=1,NSAM CRAD2 = (FI1+SWITCHX*FLOAT(K-NSAMR)**2)*HW IF (CRAD2.LT.-50.0) THEN BUF(K) = B ELSE BUF(K) = B+(BUF(K)-B)*EXP(CRAD2) ENDIF ENDDO CALL WRTLIN(LUNO,BUF,NSAM,I+(J-1)*NROW) ENDDO ENDDO ELSE C DISK MASKING DO J = 1,NSLICE FI2 = SWITCHZ*FLOAT(J-NSLICER)**2 DO I = 1,NROW FI1 = SWITCHY*FLOAT(I-NROWR)**2+FI2 CALL REDLIN(LUNI,BUF,NSAM,I+(J-1)*NROW) DO K = 1,NSAM CRAD2 = FI1 + SWITCHX*FLOAT(K-NSAMR)**2 IF (CRAD2.GT.RAD2 .OR. CRAD2.LT.RADI2) BUF(K) = B IF (ANS .EQ. 'C' .AND. CRAD2 .LT. RADI2) BUF(K) = BI ENDDO CALL WRTLIN(LUNO,BUF,NSAM,I+(J-1)*NROW) ENDDO ENDDO ENDIF END