C ++******************************************************************** C * C HKMC.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 * C * C PURPOSE: * C * C PARAMETERS: * C C IMAGE_PROCESSING_ROUTINE C * C 0 2 3 4 5 6 7 * C23456789012345678901234567890123456789012345678901234567890123456789012 C*********************************************************************** SUBROUTINE HKMC & (CIRSEED,CIRC,CIRNEW,NK,LCIRC,IP,IQ,ES,EC,E,DIST,ROT,NRING, & TEMP,MAXRIN,JACUP,NUMR,MAXIT,NIMA,MODE,TOTAL,NOUT,MIRROR) INTEGER NUMR(3,NRING),MAXRIN INTEGER*2 IQ(NK),IP(NIMA) DIMENSION CIRC(LCIRC,NIMA),CIRSEED(LCIRC,NK),CIRNEW(LCIRC) DIMENSION DIST(NIMA),ROT(NIMA) DOUBLE PRECISION TEMP(MAXRIN,2), & ENER,PEAK,PEAM,SNEW,EC(NIMA),ES(NK),E(NK) DOUBLE PRECISION F,DS,WSS,WSSMAX,TOTAL LOGICAL*1 CH_ANG,ZEROGR CHARACTER*1 MODE,MIRROR C C ES - ENERGY OF A SEED C E - SUM OF SQUARES IN A CLUSTER C DO J=1,NK ES(J)=ENER(CIRSEED(1,J),LCIRC,NRING,NUMR,MODE) ENDDO C WSSMAX=1.0D30 C DO ITER=1,MAXIT WRITE(NOUT,2020) ITER 2020 FORMAT(' Iteration #',I3) C C FIND ASSIGNMENT C IQ=0 DO I=1,NIMA F=1.0D30 DO J=1,NK CALL CROSRNG(CIRSEED(1,J),CIRC(1,I), & LCIRC,NRING,TEMP,TEMP(1,2),MAXRIN,JACUP,NUMR,PEAK,TOT,MODE) IF(MIRROR.EQ.'M') THEN CALL CROSRMG(CIRSEED(1,J),CIRC(1,I), & LCIRC,NRING,TEMP,TEMP(1,2),MAXRIN,JACUP,NUMR,PEAM,TMT,MODE) IF(PEAM.GT.PEAK) THEN DS=ES(J)+EC(I)-2.0*PEAM C IF(DS.LT.F) THEN ROT(I)=-TMT F=DS IR=J ENDIF GOTO 151 ENDIF ENDIF DS=ES(J)+EC(I)-2.0*PEAK C WRITE(NOUT,*) 'Object #',I,' Group #',J,' Dist.=',DS IF(DS.LT.F) THEN ROT(I)=TOT F=DS IR=J ENDIF 151 CONTINUE ENDDO IP(I)=IR IQ(IR)=IQ(IR)+1 C WRITE(NOUT,*) 'Set to group #',IP(I) ENDDO WRITE(NOUT,*) ' Number of objects in clusters' LK=NK/10 MK=MOD(NK,10) IP1=1 IP2=0 DO J=1,LK+MIN0(MK,1) IP2=MIN0(IP1+9,NK) WRITE(NOUT,2092) (I,I=IP1,IP2) WRITE(NOUT,2092) (IQ(I),I=IP1,IP2) WRITE(NOUT,2092) IP1=IP1+10 ENDDO 2092 FORMAT(1X,10I5) C WRITE(NOUT, *,' Assignment' C WRITE(NOUT, 2002,IP C C BUILD NEW CENTERS WITH CORRECTIONS C DO 201 J=1,NK IF(IQ(J).EQ.1) THEN DO I=1,NIMA II=I IF(IP(I).EQ.J) GOTO 468 ENDDO 468 CONTINUE CIRSEED(:,J)=CIRC(:,II) ROT(II)=0.0 DIST(II)=0.0 ES(J)=ENER(CIRSEED(1,J),LCIRC,NRING,NUMR,MODE) E(J)=0.0 GOTO 201 ENDIF CIRSEED(:,J)=0.0 C C BUILD FIRST AVERAGE C IQ(J)=0 DO 203 I=1,NIMA IF(IP(I).NE.J) GOTO 203 IQ(J)=IQ(J)+1 TOT=ROT(I) IMI=IQ(J) IF(TOT.GT.0.0) THEN CALL UPDTC(CIRSEED(1,J),CIRC(1,I), & LCIRC,NRING,NUMR,TOT,MAXRIN,IMI) ELSE CALL UPDTM(CIRSEED(1,J),CIRC(1,I), & LCIRC,NRING,NUMR,-TOT,MAXRIN,IMI) ENDIF 203 CONTINUE ES(J)=ENER(CIRSEED(1,J),LCIRC,NRING,NUMR,MODE) E(J)=0.0 C C2001 FORMAT(10(1X,F6.1)) C C ITERATIONS TO GET BETTER APPROXIMATION C ITR=0 901 CONTINUE ITR=ITR+1 CH_ANG=.FALSE. SNEW=0.0 IQT=0 C IF(MOD(ITR,2).EQ.1) THEN K1=NIMA K2=1 K3=-1 ELSE K1=1 K2=NIMA K3=1 ENDIF DO 812 IMI=K1,K2,K3 IF(MOD(ITR,2).EQ.1) THEN KTN=K1-IMI+1 ELSE KTN=IMI ENDIF IF(IP(IMI).NE.J) GOTO 812 IQT=IQT+1 CALL CROSRNG & (CIRSEED(1,J),CIRC(1,IMI),LCIRC,NRING,TEMP,TEMP(1,2), & MAXRIN,JACUP,NUMR,PEAK,TOT,MODE) IF(MIRROR.EQ.'M') THEN CALL CROSRMG & (CIRSEED(1,J),CIRC(1,IMI),LCIRC,NRING,TEMP,TEMP(1,2), & MAXRIN,JACUP,NUMR,PEAM,TMT,MODE) IF(PEAM.GT.PEAK) THEN IF(ROT(IMI) .NE. -TMT) THEN CH_ANG=.TRUE. ROT(IMI)=-TMT ENDIF F=ES(J)+EC(IMI)-2.0*PEAM CALL UPDTM(CIRNEW,CIRC(1,IMI),LCIRC,NRING,NUMR,TMT,MAXRIN,IQT) GOTO 152 ENDIF ENDIF IF(ROT(IMI).NE.TOT) THEN CH_ANG=.TRUE. ROT(IMI)=TOT ENDIF F=ES(J)+EC(IMI)-2.0*PEAK CALL UPDTC(CIRNEW,CIRC(1,IMI),LCIRC,NRING,NUMR,TOT,MAXRIN,IQT) 152 SNEW=SNEW+F DIST(IMI)=F 812 CONTINUE c WRITE(NOUT,5515) j,itr,snew,CH_ANG c5515 format(' Group #',i2,' Iteration #',i4,' SNEW=',1pe13.6,2x,l1) IF((E(J).EQ.0.0.OR.(SNEW.LE.E(J))).AND.CH_ANG) THEN CIRSEED(:,J)=CIRNEW ES(J)=ENER(CIRNEW,LCIRC,NRING,NUMR,MODE) E(J)=SNEW GOTO 901 ENDIF E(J)=SNEW 201 CONTINUE ZEROGR=.FALSE. DO I=1,NK IF(IQ(I).EQ.0) THEN ZEROGR=.TRUE. DMAX=-1.E30 DO J=1,NIMA IF(DIST(J).GT.DMAX) THEN DMAX=DIST(J) MAX=J DIST(J)=-1.E30 ENDIF ENDDO CIRSEED(:,I)=CIRC(:,MAX) WRITE(NOUT,5992) I 5992 FORMAT(' WARNING! EMPTY CLUSTER DETECTED #',I3) WRITE(NOUT,5993) MAX 5993 FORMAT(' NEW SEED CREATED FROM OBJECT #',I5) ENDIF ENDDO IF(.NOT.ZEROGR) THEN C C CALCULATE WITHIN SUM OF SQUARES C WSS=SUM(E) WNK=SUM(E/IQ) C WRITE(NOUT,2060)WSS 2060 FORMAT(' WITHIN SUM OF SQUARES=',1PE10.3) WRITE(NOUT,2008) E 2008 FORMAT(' WITHIN-GROUP SUM OF SQUARES:',/,5(2X,1PD10.3)) C WRITE(NOUT, 2001,(ANG(ROT(J),MODE),J=1,NIMA) C C HALTING RULE C IF(ABS(WSS-WSSMAX)/WSSMAX.LT.0.001) EXIT WSSMAX=WSS ENDIF C END OF LOOP OVER ITERATIONS ENDDO IF(ITER.GE.MAXIT) & WRITE(NOUT,*) ' MAXIMUM NUMBER OF ITERATIONS REACHED' WRITE(NOUT,*) ' NUMBER OF OBJECTS IN CLUSTERS' LK=NK/10 MK=MOD(NK,10) IP1=1 IP2=0 DO J=1,LK+MIN0(MK,1) IP2=MIN0(IP1+9,NK) WRITE(NOUT,2092)(I,I=IP1,IP2) WRITE(NOUT,2092)(IQ(I),I=IP1,IP2) WRITE(NOUT,2092) IP1=IP1+10 ENDDO WRITE(NOUT,*)' ASSIGNMENT' WRITE(NOUT,2002)IP 2002 FORMAT(1X,20I3) WRITE(NOUT, 2003)WSS,TOTAL,TOTAL-WSS 2003 FORMAT(' SUM OF SQUARES:',/, & ' WITHIN =',1PD10.3,' TOTAL=',1PD10.3,' BETWEEN=',1PD10.3) WRITE(NOUT,2004)WSS*(TOTAL-WSS),NK*WSS,WNK 2004 FORMAT(' CRITERIA: W*B, K*W, SUM(WK/NK)',/, & 6X,3(5X,1PE10.3)) END