C++********************************************************************* C C NOISE.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 NOISE: CALCULATE THE BACKGROUND NOISE OF POWER SPECTRUM AND SUBTRACT IT C USING LEAST-SQUARE METHOD TO FIT POINTS INTO GAUSSIAN PROFILE C USING STEEPEST DESCENT METHOD C F(XI,A)=A1EXP(-(XI/A2)**2)+A3 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE NOISE(IRTFLG) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' c INTEGER, PARAMETER :: NUMMIN = 120 REAL :: Y1(512),Y2(512) REAL :: X(120) REAL :: KFR(120),K(120) REAL :: KM,KS CHARACTER(LEN=MAXNAM) :: OUTNAME,IMFILE c COMMON /IOBUF/ BUF(NBUFSIZ) real :: BUF(5000) LUN1 = 8 c LUN2 = 10 c INUM = 1 MAXIM = 0 CALL OPFILEC(0,.TRUE.,IMFILE,LUN1,'O',IFORM,NSAM,NROW,NSLICE, & MAXIM,'IMAGE',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN WRITE(NOUT,10 )NSAM,NROW 10 FORMAT(' FILE DIMENSIONS:', I5,' X',I5) CALL DEFO003(INUM,N,KFR,Y2,NSAM,SPMAX,LUN1,BUF,NUMMIN,IRTFLG) IF (IRTFLG .NE. 0) RETURN CALL RDPRMI(NUM,NDUM,NOT_USED,'HOW MANY POINTS DO YOU WANT?') IF (NUM .EQ. 0) RETURN KM = SPMAX KS = KM / FLOAT(NSAM) WRITE(NOUT,*) ' INPUT THE POINT YOU WANT:' DO I=1,NUM WRITE(NOUT,40) I 40 FORMAT(' POINT #',I4,': ') CALL RDPRM(K(I),NOT_USED,'R [POINT]=?') C GET VALUE OF Y1, BILINEAR INTERPRELATAION Y1(I) = (1.-(K(I)-INT(K(I)))) * BUF(INT(K(I)))+ & (K(I)-INT(K(I))) * BUF(INT(K(I))+1) c WRITE(NOUT,*) K(I),Y1(I) c loc = int(k(i)) c tbuf = buf(loc) c write(NOUT,*) '----------',K(I), loc,'=', buf,':',Y1(I) ENDDO C GET X(I) VALUE DO I=1,NUM X(I) = K(I) * KS ENDDO C START HACK TO AVOID JUMP INTO LOOP --------------- C FOR SINGLE INPUT POINT IF (NUM .EQ. 1) THEN CALL RDPRM(A2,NOT_USED,'SUPPLY A2 VALUE [A-1]:') A3 = 0 A1 = Y1(1) / (EXP(-(X(1)/A2)**2)) P = A2 / KS WRITE(NOUT,*)'A1=',A1,' A2=',A2,'(A-1) ', & P,'POINTS',' A3=',A3 MAXIM = 0 CALL OPFILEC(LUN1,.TRUE.,OUTNAME,LUN2,'U',IFORM, NSAM,1,1, & MAXIM,'NOISE DELETED PROFILE',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN MAXIM = 0 CALL OPFILEC(0,.FALSE.,IMFILE,LUN1,'O',IFORM, & NSAM,NROW,NSLICE, & MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN CALL REDLIN(LUN1,Y2,NSAM,1) C ADD AN ADDITIONAL VALUE, SO ALL THE MINIMI ARE ABOVE ZERO, AT THE C OTHER PLACES VALUE BELOW ZERO ARE SET TO ZERO. XMIN = 0 DO I=1,NSAM X1 = FLOAT(I) * KS Y1(I) = A1 * EXP(-(X1/A2)**2)+A3 Y2(I) = Y2(I) - Y1(I) DO J=1,NUM IF (I .EQ. INT(K(J))-1 .OR. I .EQ. INT(K(J)) & .OR. I .EQ. INT(K(J))+1) THEN IF (XMIN .GT. Y2(I)) XMIN = Y2(I) ENDIF ENDDO ENDDO DO I=1,NSAM Y2(I) = Y2(I)-XMIN IF (Y2(I) .LT. 0) Y2(I) = 0 ENDDO CALL WRTLIN(LUN2,Y2,NSAM,1) CLOSE(LUN1) CLOSE(LUN2) RETURN ENDIF C END HACK TO AVOID JUMP THAT CAUSES COMPILE PROBLEM ----------- C SET INITIAL VALUE FOR A1,A2 AND A3 IF (NUM .GT. 2) THEN A3 = Y1(NUM)*0.8 A2 = SQRT((X(2)**2-X(1)**2)/LOG(Y1(1)/Y1(2))) A1 = (Y1(1)-A3)/EXP(-(X(1)/A2)**2) c WRITE(6,*) 'x(1), x(2) :',x(1), x(2) c WRITE(6,*) 'BUF(1), BUF(2) :',BUF(1), BUF(2) c WRITE(6,*) 'Y1(1), Y1(2) :',Y1(1), Y1(2) c write(6,*) 'Y1(1)/Y1(2) :',Y1(1)/Y1(2) c bot = LOG(Y1(1) / Y1(2)) c write(6,*) 'bot:',bot ELSEIF(NUM .GT. 1) THEN A3 = 0 A2 = SQRT((X(2)**2-X(1)**2)/LOG(Y1(1)/Y1(2))) A1 = (Y1(1)-A3) / EXP(-(X(1)/A2)**2) ENDIF c write(6,*) 'initial a1 :',a1 c write(6,*) 'initial a2 :',a2 c write(6,*) 'initial a3 :',a3 C SET INITIAL VALUE OF X**2 DO I=1,NUM Y2(I) = A1 * EXP(-(X(I)/A2)**2)+A3 ENDDO X0=0 DO I=1,NUM X0 = X0+(Y1(I)-Y2(I))**2 ENDDO NSTEP = 0 999 DA1 = 0.001 * A1 DA2 = 0.001 * A2 DA3 = 0.001 * A3 C CALCULATE dX**2 / dA1 DO I=1,NUM Y2(I) = (A1 + 0.1 * DA1) * EXP(-(X(I)/A2)**2) + A3 ENDDO X1 = 0 DO I=1,NUM X1 = X1 + (Y1(I) - Y2(I))**2 ENDDO DXA1 = (X1-X0) / (0.1*DA1) C CALCULATE DERIVATE dX**2/dA2 DO I=1,NUM Y2(I) = A1 * EXP(-(X(I)/(A2+0.1*DA2))**2)+A3 ENDDO X1 = 0 DO I=1,NUM X1 = X1+(Y1(I)-Y2(I))**2 ENDDO DXA2 = (X1-X0)/(0.1*DA2) C CALCULATE DERIVATE dX**2/dA3 DO I=1,NUM Y2(I) = A1*EXP(-(X(I)/A2)**2)+(A3+0.1*DA3) ENDDO X1 = 0 DO I=1,NUM X1 = X1+(Y1(I)-Y2(I))**2 ENDDO IF (NUM .LT. 3) THEN DXA3 = 0 ELSE DXA3 = (X1-X0)/(0.1*DA3) ENDIF C .......................... SUM = SQRT((DXA1*DA1)**2+(DXA2*DA2)**2+(DXA3*DA3)**2) A1 = A1-DXA1*DA1**2/SUM A2 = A2-DXA2*DA2**2/SUM A3 = A3-DXA3*DA3**2/SUM C CRITERIA FOR ITERATION ...... CALCULATE THE Y2 DO I=1,NUM Y2(I) = A1*EXP(-(X(I)/A2)**2)+A3 ENDDO X1 = 0 DO I=1,NUM X1 = X1+(Y1(I)-Y2(I))**2 ENDDO D = X1-X0 C SOLVE GOTO 899 BY REWRITING THE IF () THEN .. LOOP . ML 7/5/95 IF (D .LE. 0) THEN X0 = X1 P = A2 / KS C WRITE(NOUT,*)'A1=',A1,' A2=',A2,'(A-1) ',P,'POINTS',' A3',A3 NSTEP = NSTEP + 1 GOTO 999 ENDIF A1 = A1 + 0.5*DXA1*DA1**2/SUM A2 = A2 + 0.5*DXA2*DA2**2/SUM A3 = A3 + 0.5*DXA3*DA3**2/SUM P = A2 / KS WRITE(NOUT,*)'A1=',A1,' A2=',A2,'(A-1) ', & P,'POINTS',' A3=',A3 MAXIM = 0 CALL OPFILEC(LUN1,.TRUE.,OUTNAME,LUN2,'U',IFORM, NSAM,1,1, & MAXIM,'NOISE DELETED PROFILE',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN MAXIM = 0 CALL OPFILEC(0,.FALSE.,IMFILE,LUN1,'O',IFORM, & NSAM,NROW,NSLICE, & MAXIM,' ',.FALSE.,IRTFLG) IF (IRTFLG .NE. 0) RETURN CALL REDLIN(LUN1,Y2,NSAM,1) C ADD AN ADDITIONAL VALUE, SO ALL THE MINIMI ARE ABOVE ZERO, AT THE C OTHER PLACES VALUE BELOW ZERO ARE SET TO ZERO. XMIN=0 DO I=1,NSAM X1 = FLOAT(I)*KS Y1(I) = A1*EXP(-(X1/A2)**2)+A3 Y2(I) = Y2(I)-Y1(I) DO J=1,NUM IF ( I .EQ. INT(K(J))-1 .OR. I .EQ. INT(K(J)) & .OR. I .EQ. INT(K(J))+1) THEN IF (XMIN .GT. Y2(I)) XMIN = Y2(I) ENDIF ENDDO ENDDO DO I=1,NSAM Y2(I) = Y2(I)-XMIN IF (Y2(I) .LT. 0) Y2(I) = 0 ENDDO CALL WRTLIN(LUN2,Y2,NSAM,1) CLOSE(LUN1) CLOSE(LUN2) RETURN END