C++********************************************************************* C C SPEAKC.F REVISED 5/21/86 MR C RENAMED FROM SPEAKG, I2 8/18/87 MR C SPEAKC.F USED REG_GET AUG 00 ARDEAN LEITH C PARAMETERS FEB 2001 ARDEAN LEITH C NPTR MOVED AWAY FROM BUF SEP 2003 ARDEAN LEITH C IPAR BUG SEP 2007 ARDEAN LEITH 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 PURPOSE C C THIS SUBROUTINE SEARCHES FOR THE ML HIGHEST PEAKS IN THE (REAL) C IMAGE FILNAM AND PRINTS OUT POSITIONS AND VALUES OF THESE PEAKS. C HAS SUB-PIXEL ACCURACY USING BOTH CG AND POLYGONAL WEIGHTING C EXCLUSIVELY C C SPEAKC(FILNAM,LUN,NSAM,NROW,MAXDIM,OPTO,NDOC,ML,NOR) C FILNAM FILE NAME C LUN LOGICAL UNIT NUMBER OF IMAGE C NSAM,NROW DIMENSIONS OF IMAGE C MAXDIM MAXIMUM BUFFER SPACE AVAILABLE C OPTO OUTPUT OPTION C OPTO=' ' DEFAULT: NO DOCUMENT OUTPUT C OPTO='X' I.E.,REGISTER LIST FOLLOWING): C OUTPUT OF POSITION & VALUE OF PEAK IN REGISTERS C OPTO='D' DOCUMENT OUTPUT: NUMBER,POSITION, AND VALUE C OF PEAK ARE WRITTEN INTO A DOCUMENT FILE C NDOC LOGICAL UNIT NUMBER FOR DOCUMENT FILE C ML NUMBER OF PEAKS WANTED C NOR NEW ORIGIN WANTED C C REGISTER POSITIONS C 1= CG X-SHIFT C 2= CG Y-SHIFT C 3= ABSOLUTE PEAK HEIGHT C 4= RATIO PEAK HEIGHT C 5= SUB-PIXEL NON-CG FLOATING PT. X-SHIFT C 6= SUB-PIXEL NON-CG FLOATING PT. Y-SHIFT C 7= VALUE OF NON-CG EXTREMUM OF SUB-PIXEL PARABOLOID C C DOC FILE POSITIONS FOR 'PK DC' C 1= CG X-SHIFT C 2= CG Y-SHIFT C 3= ABSOLUTE PEAK HEIGHT C 4= RATIO PEAK HEIGHT C 5= NEGATIVITY WARNING FLAG C C23456789012345678901234567890123456789012345678901234567890123456789012 C--********************************************************************* SUBROUTINE SPEAKC(FILNAM,LUN,NSAM,NROW,MAXDIM,OPTO,NDOC,ML,NOR) INCLUDE 'CMBLOCK.INC' INCLUDE 'CMLIMIT.INC' COMMON /IOBUF/ BUF(NBUFSIZ) COMMON NPTR(64) INTEGER, DIMENSION(ML) :: IPOS REAL, DIMENSION(ML) :: ILIST,KLIST, KXLIST,IXLIST REAL, DIMENSION(ML) :: ALIST,RLIST REAL, DIMENSION(9) :: RSQ REAL DLIST(6) CHARACTER *(*) FILNAM CHARACTER OPTO,POSE,NULL #ifdef USE_MPI include 'mpif.h' ICOMM = MPI_COMM_WORLD MPIERR = 0 CALL MPI_COMM_RANK(ICOMM, MYPID, MPIERR) #else MYPID = -1 #endif NULL = CHAR(0) IF (MAXDIM .LT. 3*NSAM+200) THEN CALL ERRT(6,'SPEAKC',NE) RETURN ENDIF DO K = 1,ML IPOS(K) = 0 ILIST(K) = 0.0 KLIST(K) = 0.0 RLIST(K) = 0.0 ALIST(K) = -HUGE(XVAL) ENDDO NTAB = 1 NXCTR = NSAM/2+1 NYCTR = NROW/2+1 IF (NOR .NE. 0) THEN CALL RDPRMI(NXCTR,NYCTR,NOT_USED,'NEW ORIGIN COORDINATES') CALL RDPRMI(NTAB,NDUM,NOT_USED, & 'ENTER PEAK NUMBER FOR RATIO') IF (NTAB .GT. ML) THEN CALL ERRT(25,'SPEAK',NE) RETURN ENDIF ENDIF 30 CALL RDPRM2(CG,CG2,NOT_USED, & 'ELLIPSE AXES (X,Y) FOR CGR CALCULATION') IF (CG2 .EQ. 0) CG2 = CG CG2SQ = CG2*CG2 CGSQ = CG*CG CALL RDPRMC(POSE,NC,.TRUE.,'POSITIVITY ENFORCED? (Y/N)', & NULL,IRTFLG) CALL RDPRM(FNEIGH,NOT_USED,'NEIGHBORHOOD DISTANCE') CALL RDPRMI(IEWX,IEWY,NOT_USED,'EDGE EXCLUSION WIDTH X, Y') IF (IEWY .LE. 0) IEWY = IEWX C MAKE SURE THAT THE EDGE EXCLUSION IS AT LEAST AS LARGE AS THE C CENTER OF GRAVITY AXES: IF (IEWX.LT.CG) IEWX = CG IF (IEWY.LT.CG2) IEWY = CG2 WRITE(NOUT,2001) IEWX,IEWY 2001 FORMAT(' EDGE EXCLUSION WIDTH USED: X: ',I4,' Y: ',I4) IESAM = NSAM - IEWX IEROW = NROW - IEWY C CALCULATE HOW MANY RECORDS HAVE TO BE READ IN TO COVER THE COMPLETE C FIELD USED IN THE CENTER OF GRAVITY CALCULATIONS: NREC = CG2+0.5 C MAKE DIAMETER AN ODD NUMBER (=TOTAL NUMBER OF RECORDS TO C BE AVAILABLE AT ANY TIME) NREC = NREC*2+1 IF (CG .EQ. 0.0) NREC=3 NMAX = 0 NEG = 0 C INITIALIZE MREAD CALL MREAD(0,BUF,NSAM,NREC,NPTR) NREC2 = NREC/2+1 NRECH = NREC/2 DO I = NREC2,NROW-NREC2+1 CALL MREAD(LUN,BUF,NSAM,NREC,NPTR) C GET BEGINNING LOCATION OF THREE LINES I1A = (NPTR(NREC2-1)-1)*NSAM I2A = (NPTR(NREC2)-1)*NSAM I3A = (NPTR(NREC2+1)-1)*NSAM DO K = NREC2,NSAM-NREC2+1 A = BUF(K+I2A) C CHECK IF A IS SMALLER THAN THE 8 SURROUNDING POINTS: IF (A .LE. BUF(K-1+I2A)) CYCLE IF (A .LE. BUF(K-1+I1A)) CYCLE IF (A .LE. BUF(K+I1A)) CYCLE IF (A .LE. BUF(K+1+I1A)) CYCLE IF (A .LE. BUF(K+1+I2A)) CYCLE IF (A .LE. BUF(K+1+I3A)) CYCLE IF (A .LE. BUF(K+I3A)) CYCLE IF (A .LE. BUF(K-1+I3A)) CYCLE C MAKE SURE THAT PEAK NOT NEAR EDGE (DEFINED BY IEW) IF (K.LE.IEWX .OR. K.GE.IESAM) CYCLE IF (I.LE.IEWY .OR. I.GE.IEROW) CYCLE C NEW LOCAL PEAK FOUND NMAX = NMAX + 1 c write(0,*) ' nmax: ',nmax,' ml: ',ml,' A: ',A c write(0,*) ' ' C CHECK IF PEAK IS LARGER THAN THE ML PEAKS FOUND PREVIOUSLY C AND PUT IT INTO THE CORRECT RANK POSITION DO L = 1,ML !L IS PEAK NUMBER IF (A .LE. ALIST(L)) CYCLE ! CYCLE IF SMALLER IF (L .NE. ML) THEN C NOT LAST PEAK NUMBER, REORDER THE ML PEAKS L1 = L+1 DO J = L1,ML JO = ML-J+L1 JN = ML-J+L1-1 ALIST(JO) = ALIST(JN) ILIST(JO) = ILIST(JN) KLIST(JO) = KLIST(JN) IPOS(JO) = IPOS(JN) ENDDO ! END OF: DO J = L1,ML ENDIF ! END OF: IF (L .NE. ML) THEN C STORE PEAK HEIGHT AND COORDINATES IF (L .EQ. 1) THEN IPAR = I KPAR = K ENDIF ALIST(L) = A IPOS(L) = 0 ILIST(L) = I KLIST(L) = K C CENTER OF GRAVITY SEARCH IF (CG .EQ. 0.0) EXIT CNY = 0 CNX = 0 CGR = 0 PCORR = 0. IF (POSE .EQ. 'Y') THEN C POSITIVITY OF AREA IS ENFORCED, DETERMINE AREA'S MINIMUM: AMIN = 1.E10 DO ICG=1,NREC IND1 = (NPTR(ICG)-1)*NSAM YSQ = (ICG-NREC2)**2 DO KCG=K-NRECH,K+NRECH XSQ = (KCG-K)**2 CRIT = XSQ/CGSQ+YSQ/CG2SQ IF (CRIT .LE. 1) THEN IND2 = IND1+KCG XXX = BUF(IND2) IF (AMIN .GT. XXX) AMIN = XXX ENDIF ENDDO ! END OF: DO KCG=K-NRECH,K+NRECH ENDDO ! END OF: DO ICG=1,NREC PCORR = AMIN ENDIF ! END OF: IF (POSE .EQ. 'Y') THEN C CALCULATE THE CENTER OF GRAVITY: DO ICG=1,NREC IND1 = (NPTR(ICG)-1) * NSAM YSQ = (ICG-NREC2) ** 2 DO KCG=K-NRECH,K+NRECH XSQ = (KCG-K) **2 CRIT = XSQ / CGSQ + YSQ / CG2SQ IF (CRIT .GT. 1) CYCLE C POINT INSIDE ELLIPTICAL AREA FOR CGR DETERMINATION: IND2 = IND1 + KCG ADD = BUF(IND2) - PCORR IF (ADD .LT. 0) THEN IPOS(L) = 1 NEG = NEG + 1 CCCC WRITE(NOUT,1999) ADD,ICG,KCG 1999 FORMAT(' *** NEGATIVE VALUE ',G12.4, & ' FOUND AT PIXEL ',2I5,' CGR MEANINGLESS') ENDIF CGR = CGR + ADD CNY = CNY + FLOAT(ICG-NREC2)*ADD CNX = CNX + FLOAT(KCG-K)*ADD ENDDO ! END OF: DO KCG=K-NRECH,K+NRECH ENDDO ! END OF: DO ICG=1,NREC IF (IPOS(L) .EQ. 0) THEN C MEANINGFULL CGR ILIST(L) = CNY / CGR + ILIST(L) KLIST(L) = CNX / CGR + KLIST(L) EXIT ENDIF ENDDO ! END OF: DO L = 1,ML ENDDO ! END OF: DO K = NREC2,NSAM-NREC2+1 ENDDO IF (NMAX .EQ. 0) THEN IF (MYPID .LE. 0) WRITE(NDAT,*) '*** NO PEAKS FOUND' IF (MYPID .LE. 0 .AND. NOUT .NE. NDAT) & WRITE(NOUT,*) '*** NO PEAKS FOUND' CALL REG_SET_NSEL(1,5,0.0, 0.0, 0.0, 0.0, 0.0,IRTFLG) CALL REG_SET_NSEL(6,2,0.0, 0.0, 0.0, 0.0, 0.0,IRTFLG) RETURN ENDIF IF (MYPID .LE. 0) THEN IF (NEG .GT. 0) THEN WRITE(NOUT,304)NEG 304 FORMAT(/,' WARNING: ',I6,' NEGATIVE PIXEL VALUES FOUND', & ' DURING CALCULATION OF CGR',/) ENDIF WRITE(NDAT,299) IF (NDAT .NE. NOUT) WRITE(NOUT,299) 299 FORMAT(' NO NSAM NROW X Y VALUE', & ' RATIO',5X,' 1 IF NEGATIVE'/5X,58X,'VALUES IN CGR AREA') ENDIF MLIST = MIN(ML,NMAX) C FOR EACH PEAK, MAKE SURE IT DOES NOT FALL WITHIN RADIUS=FNEIGH C OF A HIGHER-RANKING ONE FN2 = FNEIGH * FNEIGH C CHANGED 2/8/88 MR DO L=MLIST,2,-1 FK = KLIST(L) FI = ILIST(L) DO LI=1,L-1 C IF OUTIDE EXCLUSION ZONE, CYCLE IF ((FK-KLIST(LI))**2+(FI-ILIST(LI))**2 .GT. FN2) CYCLE ALIST(L) = -HUGE(XVAL) ENDDO ENDDO NPKCNT = 0 C NOW SELECT PEAKS THAT COMPLY WITH CONDITIONS AND ENTER IN DOC-FILE DO L=1,MLIST IF (ALIST(L) .EQ. -HUGE(XVAL)) THEN C PEAK HAS BEEN EXCLUDED IF (L .NE. NTAB) CYCLE NTAB = NTAB + 1 IF (MYPID .LE. 0) WRITE(NOUT,319) NTAB IF (MYPID .LE. 0 .AND. NOUT .NE. NDAT) & WRITE(NOUT,319) NTAB 319 FORMAT(' NUMBER OF REF. PEAKS INCREASED BY ONE TO: ',I5) CYCLE ENDIF NPKCNT = NPKCNT + 1 KXLIST(L) = KLIST(L) - NXCTR IXLIST(L) = ILIST(L) - NYCTR RLIST(L) = ALIST(L) / ALIST(NTAB) IF (MYPID .LE. 0) THEN C WRITE TO TERMINAL AND RESULTS FILE WRITE(NDAT,301) NPKCNT,KLIST(L),ILIST(L), & KXLIST(L),IXLIST(L),ALIST(L),RLIST(L),IPOS(L) 301 FORMAT(I6, 4F8.2, G12.3, F10.5,5X,I2) IF (NOUT .NE. NDAT) THEN WRITE(NOUT,301) NPKCNT,KLIST(L),ILIST(L), & KXLIST(L),IXLIST(L), ALIST(L),RLIST(L),IPOS(L) ENDIF ENDIF ! END OF: IF (MYPID .LE. 0) THEN IF (OPTO .EQ. 'D') THEN C WRITE PEAK TO DOCUMENT FILE DLIST(1) = L DLIST(2) = KXLIST(L) DLIST(3) = IXLIST(L) DLIST(4) = ALIST(L) DLIST(5) = RLIST(L) DLIST(6) = IPOS(L) CALL SAVD(NDOC,DLIST,6,IRTFLG) ENDIF ENDDO ! END OF: DO L=1,MLIST CALL SAVDC CLOSE(NDOC) C SET RESULTS IN COMMAND LINE REGISTERS 1..4 CALL REG_SET_NSEL(1,4,KXLIST(NTAB),IXLIST(NTAB), & ALIST(NTAB), RLIST(NTAB),0.0,IRTFLG) C 9/25/81 PARABOLIC FIT TO THE 3X3 NEIGHBORHOOD OF THE PEAK POINT C PROGRAM SECTION SENT BY M.VAN HEEL, MODIFIED FOR SPIDER USE. JF C APPLIED ONLY TO HIGHEST RANKING PEAK KL = KPAR DO I=1,3 if (Ipar .gt. nrow .or. ipar .le. 0) then write(6,*) ' IPAR: ',IPAR,' > NROW:',nrow write(6,*) ' I: ',I CALL ERRT(101,'BAD ROW REFERENCED',NE) endif IROW = IPAR + I - 2 ! GET ROW NUMBER IF (IROW .LT. 1) IROW = IROW + NROW IF (IROW .GT. NROW) IROW = IROW - NROW C READ ORIGINAL DATA CALL REDLIN(LUN,BUF,NSAM,IROW) I1 = (I-1)*3 DO K=1,3 ISAM = KL+K-2 IF (ISAM .LT.1) ISAM=ISAM+NSAM IF (ISAM .GT. NSAM) ISAM=ISAM-NSAM C VALUES FOR NINE LOCATIONS AROUND THE PEAK & PEAK RSQ(I1+K)= BUF(ISAM) ENDDO ENDDO C FIND XSH & YSH (RANGE -1.0 .... 1.0) CALL PARABL(RSQ,XSH,YSH,PEAKV) C SET SUB-PIXEL RESULTS IN COMMAND LINE REGISTERS 5..7 CALL REG_SET_NSEL(5,3,XSH+KPAR-NXCTR, YSH+IPAR-NYCTR, & PEAKV, 0,0,IRTFLG) IF (MYPID .LE. 0) THEN WRITE(NDAT,*) ' ' WRITE(NDAT,*)' SUB PIXEL OFFSETS FOR LARGEST PEAK: ' IF (NDAT .NE. NOUT) WRITE(NOUT,*) ' ' IF (NDAT .NE. NOUT) & WRITE(NOUT,*)' SUB PIXEL OFFSETS FOR LARGEST PEAK: ' WRITE(NDAT,351)XSH,YSH,PEAKV IF (NOUT .NE. NDAT) WRITE(NOUT,351)XSH,YSH,PEAKV 351 FORMAT(' + (',F6.2,', ',F6.2,') PEAK VALUE = ',G12.5,/) IDIFF = ML - NPKCNT WRITE(NOUT,401) ML,NPKCNT,IDIFF 401 FORMAT(' **',I5,' PEAKS SPECIFIED ',I6,' PEAKS PASSED ', & I6,' PEAKS EXCLUDED (NEIGB. DIST. EDGE EXCL.)') ENDIF END