C ********************************************************************** C C ADDS.F USED FILELIST JULY 99 ARDEAN LEITH C USED F90 JULY 99 ARDEAN LEITH C SQRT(NEG) TRAP APR 00 ARDEAN LEITH C ALLOC TRAP APR 03 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 PURPOSE: ADD A SERIES OF IMAGES C C--********************************************************************* SUBROUTINE ADDS(LUN1,LUN2,LUNS,NDOC,MAXMEM) COMMON NUMB(1) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' CHARACTER (LEN=MAXNAM) :: FILNAM,FILA,FILV,FILPAT DOUBLE PRECISION :: VARAV,AV2,VV CHARACTER (LEN=1) :: NULL, SER REAL, ALLOCATABLE, DIMENSION(:) :: AVGARAY,VARARAY,INPARAY #ifndef SP_32 INTEGER *8 NVOX,IOK8 #else INTEGER *4 NVOX,IOK8 #endif NULL = CHAR(0) NUMTOT = MAXMEM -100 CALL FILELIST(.TRUE.,NDOC,FILPAT,NLETP,NUMB,NUMTOT,NUMT, & 'INPUT FILE TEMPLATE (E.G. PIC****)',IRTFLG) IF (IRTFLG .NE. 0) RETURN C GET FIRST PICTURE TO DETERMINE DIMS CALL FILGET(FILPAT,FILNAM,NLETP,NUMB(1),IRTFLG) IF (IRTFLG .NE. 0) RETURN MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILNAM,LUN2,'O',IFORM,NSAM,NROW,NSLICE, & MAXIM,' ',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN MAXIM = 0 CLOSE(LUN2) NVOX = NSAM * NROW * NSLICE C COMPLAIN IF EXCESSIVE ALLOCATION NVOX3 = NVOX * 3 CALL BIGALLOC(NVOX,IOK8,.FALSE.,.TRUE.,IRTFLG) ALLOCATE(AVGARAY(NVOX),VARARAY(NVOX),INPARAY(NVOX),STAT=IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(102,'FAILED TO ALLOCATE',NVOX3) GOTO 999 ENDIF CALL RDPRMC(SER,NUMC,.TRUE., 'ALL, ODD-EVEN (A/O)',NULL,IRT) IF (IRTFLG .NE. 0) GOTO 999 IF (SER.NE.'O' .AND. SER.NE.'E') THEN SER = 'A' LUP = 1 ELSE LUP = 2 ENDIF DO ILUP=1,LUP IF (LUP.EQ.2) THEN IF (ILUP .EQ. 1) THEN WRITE(NOUT,*) 'FOR ODD-NUMBERED IMAGES' SER = 'O' ELSE WRITE(NOUT,*) 'FOR EVEN-NUMBERED IMAGES' SER = 'E' ENDIF ENDIF CALL FILERD(FILA,NLETA,NULL,'AVERAGE',IRTFLG) IF (IRTFLG .NE. 0) GOTO 997 C COMMAND 'AS AD' MEANS THAT AVERAGE FILE ALREADY EXISTS IF (FCHAR(4:5) .EQ. 'AD') THEN C FILES ALREADY EXIST,FIND OUT HOW MANY PICTURES ALREADY AVERAGED C5 IF (NSEL(1) .EQ. 0) THEN 5 CALL REG_GET_USED(NSEL_USED) IF (NSEL_USED .LE. 0) THEN WRITE(NOUT,*) '*** MISSING OFFSET VALUE' CALL ERRT(100,'ADDS',NE) GOTO 997 ENDIF CALL RDPRMI(NUMPR,NDUM,NOT_USED, & 'ENTER NO. OF IMAGES ALREADY AVERAGED') C OPEN EXISTING AVERAGE FILE CALL OPFILEC(0,.FALSE.,FILA,LUN1,'Z',IFORM,NSAM,NROW,NSLICE, & MAXIM,' ',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 997 C OPEN OLD VARIANCE FILE CALL OPFILEC(0,.TRUE.,FILV,LUNS,'Z',IFORM,NSAM,NROW,NSLICE, & MAXIM,'VARIANCE',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 997 C GET OLD OFFSET C OFFOLD = PARAM(NSEL(1)) CALL REG_GET_NSEL(1,OFFOLD,FDUM,FDUM,FDUM,FDUM,IRTFLG) DO I=1,NROW*NSLICE CALL REDLIN(LUN1,AVGARAY((I-1)*NSAM+1),NSAM,I) ENDDO DO I=1,NROW*NSLICE CALL REDLIN(LUNS,VARARAY((I-1)*NSAM+1),NSAM,I) ENDDO VARARAY = VARARAY * (NUMPR-1) + AVGARAY * AVGARAY * NUMPR AVGARAY = AVGARAY * NUMPR ELSE C FILE DOES NOT EXIST - OPEN WITH DIMS OF FIRST FILE & INITIALIZE C OPEN NEW AVERAGE OUTPUT FILE MAXIM = 0 CALL OPFILEC(0,.FALSE.,FILA,LUN1,'U',IFORM,NSAM,NROW,NSLICE, & MAXIM,' ',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 997 C OPEN NEW VARIANCE OUTPUT FILE MAXIM = 0 CALL OPFILEC(0,.TRUE.,FILV,LUNS,'U',IFORM,NSAM,NROW,NSLICE, & MAXIM,'VARIANCE',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) GOTO 997 C CLEAR BOTH AVERAGE AND VARIANCE ARRAYS AVGARAY = 0.0 VARARAY = 0.0 NUMPR = 0 C NUMBER OF PREVIOUS PICTURES IS 0 ENDIF NUMC = 1 AV1 = 0.0 C CHANGED TO ALLOW GAPS IN FILE SERIES 7/21/89 MR.: NUMA = 0 NMT = 0 DO IFIL=1,NUMT CALL FILGET(FILPAT,FILNAM,NLETP,NUMB(IFIL),IRTFLG) IF (IRTFLG .NE. 0) THEN CALL ERRT(3,'ADDS',NE) GOTO 997 ENDIF CALL OPFILEC(0,.FALSE.,FILNAM,LUN2,'Z',IFORM, & NSAM,NROW,NSLICE, & MAXIM,' ',.TRUE.,IRTFLG) IF (IRTFLG .NE. 0) THEN C TO ALLOW GAPS IN FILE SERIES 7/21/89 MR.: WRITE(NOUT,100) FILNAM 100 FORMAT(' FILE: ',A,' NOT FOUND, FILE SKIPPED') GOTO 550 ENDIF MAXIM = 0 NMT = NMT+1 IF (SER.EQ.'O' .AND. MOD(NMT,2) .EQ. 0) GOTO 707 IF (SER.EQ.'E' .AND. MOD(NMT,2) .EQ. 1) GOTO 707 NUMA = NUMA+1 DO I=1,NROW*NSLICE CALL REDLIN(LUN2,INPARAY((I-1)*NSAM+1),NSAM,I) CONTINUE ENDDO IF (FCHAR(4:4) .NE. 'R') THEN AV2=SUM(INPARAY)/ NSAM / NROW / NSLICE INPARAY = INPARAY - AV2 ELSE AV2=0.0 ENDIF AVGARAY= AVGARAY+ INPARAY VARARAY= VARARAY+ INPARAY*INPARAY AV1 = AV1 + AV2 707 CLOSE(LUN2) 550 CONTINUE ENDDO FNUMT = NUMA + NUMPR C NOW COMPUTE NORMALIZED AVERAGE AND VARIANCE IMAGE AVGARAY = AVGARAY / FNUMT VV = DOT_PRODUCT(AVGARAY,AVGARAY) VARARAY = (VARARAY - AVGARAY*AVGARAY*FNUMT) /(FNUMT-1.0) VARAV = SUM(VARARAY) DO I=1,NROW*NSLICE CALL WRTLIN(LUN1,AVGARAY((I-1)*NSAM+1),NSAM,I) ENDDO DO I=1,NROW*NSLICE CALL WRTLIN(LUNS,VARARAY((I-1)*NSAM+1),NSAM,I) ENDDO VV = VV/REAL(NSAM*NROW*NSLICE) IF (VARAV .LT. 0.0) THEN WRITE(NOUT,*) 'WARNING: NEGATIVE TOTAL VARIANCE: ',VARAV, & ' SD SET TO ZERO!!!' IF (NOUT .NE. NDAT) & WRITE(NDAT,*) 'WARNING: NEGATIVE TOTAL VARIANCE: ',VARAV, & ' SD SET TO ZERO!!!' SAVT = 0.0 ELSE SAVT = SQRT(VARAV) ENDIF VAV = VARAV/(FLOAT(NSAM)*FLOAT(NROW)*FLOAT(NSLICE)) IF (VAV .LT. 0.0) THEN WRITE(NOUT,*) 'WARNING: NEGATIVE VARIANCE P.P. : ',VAV, & ' SD SET TO ZERO!!!' IF (NOUT .NE. NDAT) & WRITE(NDAT,*) 'WARNING: NEGATIVE TOTAL VARIANCE: ',VAV, & ' SD SET TO ZERO!!!' SAV = 0.0 ELSE SAV = SQRT(VAV) ENDIF OFF = AV1 / FLOAT(NUMA) C IF (NUMPR.NE.0.AND.NSEL(1).NE.0.AND.PARAM(NSEL(1)).NE.0.0) C & OFF = (OFFOLD*FLOAT(NUMPR)+AV1)/FLOAT(NUMA+NUMPR) CALL REG_GET_USED(NSEL_USED) IF (NUMPR.NE.0 .AND. NSEL_USED.NE.0) THEN C TRAP FOR NaN IF (OFFOLD.NE.0.0) & OFF = (OFFOLD*FLOAT(NUMPR)+AV1)/FLOAT(NUMA+NUMPR) ENDIF C IF (NSEL(1) .NE .0) PARAM(NSEL(1)) = OFF CALL REG_SET_NSEL(1,1, OFF,0.0, 0.0, 0.0, 0.0,IRTFLG) IF (NOUT .NE. NDAT) & WRITE(NDAT,7001) NUMA+NUMPR,NSAM*NROW*NSLICE,VARAV,SAVT, & VAV,SAV,OFF,VV WRITE(NOUT,7001) NUMA+NUMPR,NSAM*NROW*NSLICE,VARAV,SAVT, & VAV,SAV,OFF,VV 7001 FORMAT(' ** VARIANCE COMPUTATION BASED ON ',I7,' IMAGES ', & /' CONTAINING ',I7,' ELEMENTS', & /' TOTAL VARIANCE = ',1PG12.5,' TOTAL S.D. = ',1PG12.5, & /' VARIANCE P.P. = ',1PG12.5,' S.D. P.P. = ',1PG12.5, & /' AVERAGE OFFSET = ',1PG12.5, & /' VARIANCE OF AVERAGE IMAGE = ',1PG12.5) ENDDO 997 CLOSE(LUNS) CLOSE(LUN1) 999 IF (ALLOCATED(INPARAY)) DEALLOCATE(INPARAY) IF (ALLOCATED(VARARAY)) DEALLOCATE(VARARAY) IF (ALLOCATED(AVGARAY)) DEALLOCATE(AVGARAY) RETURN END