CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C PROGRAM IMSEGP.DET.AUTO C C VERSION 1.6 17-SEP-82 C C C C C C PROGRAM IMSEGP.DET.AUTO (IMPAC:SEGPDETA) WAS DEVELOPED C C from Program IMSEGP.DET. C C THE "IMSEGP" PROGRAMS ARE FOR SEGMENTING MULTI-CHANNEL C DIGITAL IMAGES. (FOR SINGLE-CHANNEL DATA THE "IMSEG1" PROGRAMS C MAY BE USED.) C IMSEGP.COMMON USES DISTANCE IN THE METRIC OF THE ESTIMATED C COMMON COVARIANCE MATRIX. IMSEG.DET USES DIFFERENT COVARIANCE C MATRICES, WITH ADJUSTMENT BY THE DETERMINANTS, I.E., IT USES THE C C the estimated log likelihood for the Gaussian model with C C different covariance matrices. C C C C AUTOMATIC MODE: TRIES A RANGE OF NUMBERS OF CLASSES, WITH C C AUTOMATIC SETTING OF INITIAL MEANS, EQUALLY SPACED THROUGH THE C C RANGE OF EACH VARIABLE. C C C C C C PROGRAMMED BY: C C DR. STANLEY L. SCLOVE 312/996-2681 C C DEPARTMENT OF QUANTITATIVE METHODS 312/996-2676 C C COLLEGE OF BUSINESS ADMINISTRATION C C UNIVERSITY OF ILLINOIS AT CHICAGO C C BOX 4348, CHICAGO, IL 60680 C C C C RESEARCH SUPPORTED IN PART BY: C C C C ONR CONTRACT N00014-80-C-0408, TASK NR042-443 C C ARO CONTRACT DAAG29-82-K-0155 C C C C RESTRICTIONS (CAN BE MODIFIED): C C NR, NUMBER OF ROWS OF ARRAY, AT MOST 256; C C NC, NUMBER OF COLUMNS OF ARRAY, AT MOST 256; C C IP, NUMBER OF CHANNELS, AT MOST 20; C C K, NUMBER OF CLASSES, AT MOST 29; C C ITER, MAXIMUM NUMBER OF ITERATIONS, 20. C C C C SUBROUTINE(S) CALLED: C C MATEQ, WHICH CALLS MATDT C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C DIMENSION X( 56, 56,10),SUM(29, 4) DIMENSION D(29),ICLUS( 56, 56) C DIMENSION IOTA( 56),JOTA(29) C DIMENSION TITLE(18) DIMENSION NG(29),XMEAN(29, 4) DIMENSION FMT(18) DIMENSION SS(29, 4, 4),SSD(29, 4, 4) DIMENSION WGSS( 4, 4) DIMENSION VARHAT( 4, 4),WGMS( 4, 4) DIMENSION ICLSOL( 56, 56) DIMENSION XMIN( 4),XMAX( 4) C CHANGE NO. OF CHANNELS FROM 4 TO 20 C AND NO. OF ROWS AND COLUMNS FROM 56 TO 256 C IN DIMENSION STATEMENTS WHEN ALLOWABLE REGION IS INCREASED. DIMENSION IV(20,20) DIMENSION P(20,20) C DIMENSION A(20,20) C DIMENSION ET(29) DIMENSION PG(29,20,20) C DIMENSION NT(29,29,29),IRSUM(29,29),TP(29,29,29) DIMENSION PROB(29) C C DOUBLE PRECISION SS,SUM DOUBLE PRECISION WGSS,SSD DOUBLE PRECISION VARHAT DOUBLE PRECISION P DOUBLE PRECISION DET DOUBLE PRECISION D DOUBLE PRECISION XMEAN DOUBLE PRECISION TEMPIV,TEMPJV DOUBLE PRECISION F DOUBLE PRECISION CF C DOUBLE PRECISION A C DOUBLE PRECISION ET DOUBLE PRECISION PG DOUBLE PRECISION TP,PROB C C IV IS A WORK ARRAY FOR SUBROUTINE MATEQ. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C C CONTROL CARDS: C C C C DATASET TITLE C C NUMBER OF ROWS, NR, IN FORMAT(3X,I3) C C NUMBER OF COLUMNS, NC, IN FORMAT(3X,I3) C C NUMBER OF CHANNELS, IP, IN FORMAT (3X,I2) C C FMT, IN FORMAT (18A4), E.G., (1X,F4.1) C C "FMT" WILL ALSO BE USED FOR OUTPUT: ALLOW AT LEAST ONE BLANK C C AT THE BEGINNING FOR CARRIAGE CONTROL. C C DATA, ONE NUMBER AT A TIME, IN FORMAT SPECIFIED BY FMT C C DATA IS INDEXED BY CHANNEL, ROW, AND COLUMN. C C COLUMN CHANGES FIRST, THEN ROW, THEN CHANNEL. C C C C K, NUMBER OF CLASSES, IN FORMAT (2X,I2) C C K INITIAL MEANS, IN FORMAT SPECIFIED BY FMT C C INITIAL MEANS ARE INDEXED BY CLASS AND CHANNEL. C C CHANNEL CHANGES FIRST, THEN CLASS. C C C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C READ (5,28000) TITLE C C WRITE PROGRAM INFORMATION. WRITE (6,17000) WRITE (6,53000) WRITE (6,35000) TITLE C C READ NUMBERS OF ROWS AND COLUMNS, NR AND NC. READ (5,16000) NR READ (5,16000) NC WRITE (6,33000) NR WRITE (6,34000) NC C READ NUMBER OF CHANNELS (I.E., VARIABLES), IP. READ (5,29000) IP WRITE (6,14000) IP C C READ DATA FORMAT. READ (5,28000) FMT C C READ DATA. DO 300 ICHAN = 1,IP DO 300 I = 1,NR DO 300 J = 1,NC READ (5,FMT) X(I,J,ICHAN) IF (I .EQ. 1 .AND. J .EQ. 1) GO TO 100 GO TO 200 100 CONTINUE DO 200 IVAR = 1,IP XMAX(IVAR) = X(1,1,IVAR) XMIN(IVAR) = X(1,1,IVAR) 200 CONTINUE DO 300 IVAR = 1,IP IF (X(I,J,IVAR) .LT. XMIN(IVAR)) XMIN(IVAR)= X X(I,J,IVAR) IF (X(I,J,IVAR) .GT. XMAX(IVAR)) XMAX(IVAR)= X X(I,J,IVAR) 300 CONTINUE WRITE (6,54000) WRITE (6,FMT) (XMIN(IVAR),IVAR=1,IP) WRITE (6,55000) WRITE (6,FMT) (XMAX(IVAR),IVAR=1,IP) C C READ KL AND KU, THE MIN AND MAX NUMBER OF CLASSES TO BE TRIED. READ (5,15000) KL READ (5,15000) KU WRITE (6,18000) KL,KU C C SET CONSTANTS. C PI = 3.1415927 C N = NR*NC XN = N DO 400 INTEG=1,NC IOTA(INTEG) = INTEG 400 CONTINUE C DO 6500 K = KL,KU C C C COMPUTE INITIAL MEANS: DO 500 IG=1,K DO 500 IVAR = 1,IP XMEAN(IG,IVAR) = XMIN(IVAR) + IG*(XMAX(IVAR)- X XMIN(IVAR))/(K+1) 500 CONTINUE C WRITE (6,19000) K C WRITE INITIAL MEANS C DO 600 IG = 1,K WRITE (6,20000) IG,(XMEAN(IG,IVAR),IVAR=1,IP) 600 CONTINUE C C SET CONSTANTS FOR THIS VALUE OF K. C DO 700 J = 1,K JOTA(J) = J 700 CONTINUE C C PARAMETERS FOR MODEL WITH COMMON COVARIANCE MATRIX ARE C K MEAN VECTORS AND ONE COVARIANCE MATRIX. NOPARM = K*IP + IP*(IP+1)/2 C PARAMETERS FOR MODEL WITH DIFFERENT COVARIANCE MATRICES ARE C K MEAN VECTORS AND K COVARIANCE MATRICES. NPRMDF = K*IP + K*IP*(IP+1)/2 C C NON-DISTRIBUTIONAL PARAMETERS: C K**2-BY-K TRANSITION PROB. MATRIX C GIVES (K**2)(K-1) FREE TRANSITION PROBABILITIES C NOPARM = NOPARM + (K**2)*(K-1) NPRMDF = NPRMDF + (K**2)*(K-1) C C C FOR FIRST ITERATION, MARGINAL DISTRIBUTION OF LABELS, PROB, C IS TAKEN TO BE UNIFORM. C DO 800 IG=1,K 800 PROB(IG) = 1.0/K C C ITER = 1 C ITERATIONS BEGIN HERE C 900 CONTINUE IF (ITER .EQ. 1) GO TO 1100 DO 1000 I = 1,NR DO 1000 J = 1,NC ICLSOL(I,J) = ICLUS(I,J) 1000 CONTINUE C SAVE PREVIOUS VALUE OF -2 LOG MAX LIKELIHOOD: XMN20L = XMN2LL 1100 CONTINUE C C COMMENCE DISTANCE COMPUTATIONS. C DO 2400 I = 1,NR DO 2400 J = 1,NC C INITIALIZE DISTANCES (TO BE ACCUMULATED) AT ZERO. DO 1200 L = 1,K D(L) = 0.0 1200 CONTINUE C FOR FIRST ITERATION, EUCLIDEAN DISTANCE IS USED BECAUSE C NO COVARIANCE MATRIX IS YET AVAILABLE. AFTER THE FIRST C ITERATION, DISTANCE WILL BE TAKEN IN THE METRIC OF THE C GROUP COVARIANCE MATRIX. C IF (ITER .GT. 1) GO TO 1400 DO 1300 L=1,K DO 1300 IVAR=1,IP D(L) = D(L) + ( XMEAN(L,IVAR) - X(I,J, X IVAR) )**2 1300 CONTINUE GO TO 1700 1400 CONTINUE DO 1600 L=1,K DO 1500 IVAR=1,IP TEMPIV = XMEAN(L,IVAR) - X(I,J,IVAR) DO 1500 JV=1,IP TEMPJV = XMEAN(L,JV) - X(I,J,JV) C D(L) = D(L) + TEMPIV*PG(L,IVAR,JV) X *TEMPJV C 1500 CONTINUE C 1600 CONTINUE 1700 CONTINUE C DO 2200 L = 1,K C FOR FIRST ITERATION, CLASSIFICATION IS SIMPLY BY C MINIMUM DISTANCE. OTHERWISE, THE TRANSITION C PROBABILITIES AND DETERMINANTS ENTER: IF (ITER .EQ. 1) GO TO 2200 C UP TO NOW, D(L) IS (SQUARED) DISTANCE C IT MUST NOW BE MODIFIED TO PROBABILITY, C FOR MULTIPLICATION BY THE TRANSITION PROBABILITIES. C ADD DETERMINANT OF COVARIANCE MATRIX TO DISTANCE: D(L) = D(L) + ET(L) ARG = -D(L)/2.0 IF (ARG .LT. -180.2) GO TO 1800 GO TO 1900 1800 D(L) = 0.0 GO TO 2000 1900 CONTINUE C C MOVE FROM LOG PROB SCALE TO PROB SCALE: D(L) = EXP(ARG) C D(L) IS NO LONGER A DISTANCE: C IT IS (PROPORTIONAL TO) THE BELONGING PROBABILITY 2000 CONTINUE C IF (I .EQ. 1 .OR. J .EQ. 1) GO TO 2100 IG1 = ICLUS(I,J-1) IG2 = ICLUS(I-1,J) D(L) = -TP(IG1,IG2,L)*D(L) GO TO 2200 2100 CONTINUE C CLASSIFY BOUNDARY OBSERVATIONS D(L) = -D(L)*PROB(L) 2200 CONTINUE C F = D(1) ICLUS(I,J) = 1 DO 2400 L = 2,K IF ( D(L) - F ) 2300,2400,2400 2300 F = D(L) ICLUS(I,J) = L 2400 CONTINUE WRITE (6,10000) C WRITE (6,11000) (IOTA(J), J=1,NC) DO 2500 I = 1,NR WRITE (6,12000) I, (ICLUS(I,J), J=1,NC) 2500 CONTINUE DO 2600 IG = 1,K NG(IG) = 0 DO 2600 IVAR = 1,IP SUM(IG,IVAR) = 0.0 DO 2600 JV = 1,IP SS(IG,IVAR,JV) = 0.0 SSD(IG,IVAR,JV) = 0.0 2600 CONTINUE DO 2700 I = 1,NR DO 2700 J = 1,NC IGROUP = ICLUS(I,J) NG(IGROUP) = NG(IGROUP) + 1 DO 2700 IVAR = 1,IP SUM(IGROUP,IVAR) = SUM(IGROUP,IVAR) + X(I,J, X IVAR) DO 2700 JV = 1,IP SS(IGROUP,IVAR,JV) = SS(IGROUP,IVAR,JV) + X X(I,J,IVAR)*X(I,J,JV) 2700 CONTINUE C WRITE (6,27000) (NG(IG),IG=1,K) C DO 2800 IVAR = 1,IP DO 2800 JV = 1,IP WGSS(IVAR,JV) = 0.0 2800 CONTINUE C C DO 3100 IG = 1,K IF (NG(IG) .EQ. 0) GO TO 2900 GO TO 3000 2900 WRITE (6,60000) IG GO TO 6600 3000 CONTINUE DO 3100 IVAR = 1,IP C COMPUTE MEAN VECTORS XMEAN(IG,IVAR) = SUM(IG,IVAR)/NG(IG) C COMPUTE SUM-OF-PRODUCTS MATRICES DO 3100 JV = 1,IP CF = SUM(IG,IVAR)*SUM(IG,JV)/NG(IG) SSD(IG,IVAR,JV) = SS(IG,IVAR,JV) - CF 3100 CONTINUE C C C C POOL: C DO 3200 IG = 1,K DO 3200 IVAR = 1,IP DO 3200 JV = 1,IP WGSS(IVAR,JV) = WGSS(IVAR,JV) + SSD(IG,IVAR, X JV) 3200 CONTINUE C C C COMPUTE VARHAT, MLE OF COMMON COVARIANCE MATRIX: DO 3300 IVAR = 1,IP DO 3300 JV = 1,IP VARHAT(IVAR,JV) = WGSS(IVAR,JV)/N 3300 CONTINUE C COMPUTE DET(VARHAT) IDET = 1 NRS1 = 0 DO 3400 IVAR=1,IP DO 3400 JV = 1,IP A(IVAR,JV) = VARHAT(IVAR,JV) 3400 CONTINUE CALL MATEQ(A,IP,20,JFLG,DET,IDET,IV,NRS1,P,20) C C GENERAL FORM OF CALL IS: C CALL MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C SEE SUBROUTINE LISTING FOR FULLER EXPLANATION. C DET(VARHAT) = DET*10.0**IDET C (ACTUAL DET. = DET*10**IDET) C IF (JFLG .GT. 0) WRITE (6,43000) JFLG C XIDET = IDET XLGDET = DLOG(DET) + XIDET*ALOG(10.0) XMN2LL = N*(IP*ALOG(2.0*PI) + IP + XLGDET) C C IF (ITER .EQ. 1) GO TO 3600 C C COMPARE NEW SEGMENTATION WITH OLD C DO 3500 I = 1,NR DO 3500 J = 1,NC IF (ICLUS(I,J) .EQ. ICLSOL(I,J)) GO TO 3500 GO TO 3600 3500 CONTINUE GO TO 6000 3600 CONTINUE C C IF NEW SEGMENTATION DIFFERS FROM OLD, C COMPUTE AND WRITE NEW STATISTICS C DO 3700 IG = 1,K WRITE (6,26000) IG, (XMEAN(IG,IVAR),IVAR=1,IP) 3700 CONTINUE C C COMPUTE TRANSITION PROBABILITY MATRIX: C DO 3800 I1 = 1,K DO 3800 I2 = 1,K DO 3800 J = 1,K NT(I1,I2,J) = 0 3800 CONTINUE DO 3900 I = 2,NR DO 3900 J = 2,NC IM1 = I-1 JM1 = J-1 INORTH = ICLUS(IM1,J) IWEST = ICLUS(I,JM1) IY = ICLUS(I,J) NT(IWEST,INORTH,IY) = NT(IWEST,INORTH,IY) + 1 3900 CONTINUE DO 4000 I1 = 1,K DO 4000 I2 = 1,K IRSUM(I1,I2) = 0 4000 CONTINUE DO 4100 I1 = 1,K DO 4100 I2 = 1,K DO 4100 J = 1,K IRSUM(I1,I2) = IRSUM(I1,I2) + NT(I1,I2,J) 4100 CONTINUE DO 4400 I1=1,K DO 4400 I2=1,K XDENOM=IRSUM(I1,I2) IF (XDENOM .EQ. 0.0) GO TO 4200 GO TO 4300 4200 XDENOM = K 4300 CONTINUE C DO 4400 J = 1,K XNUM = NT(I1,I2,J) C IF THERE ARE NO TRANSITIONS FROM (I1,I2), THEN TP(I1,I2,J) C IS SET EQUAL TO ZERO, FOR ALL J = 1,2,...,K. C 4400 TP(I1,I2,J) = XNUM/XDENOM 4500 CONTINUE C C COMPUTE MARGINAL DISTRIBUTION OF LABELS: DO 4600 IG = 1,K PROB(IG) = NG(IG) PROB(IG) = PROB(IG)/N 4600 CONTINUE C TRANS = 0.0 DO 4700 I1=1,K DO 4700 I2=1,K DO 4700 J=1,K ITEST = NT(I1,I2,J) IF (ITEST .EQ. 0) GO TO 4700 IF (TP(I1,I2,J) .EQ. 0.) GO TO 4700 TRANS = TRANS + NT(I1,I2,J)*DLOG(TP(I1,I2,J)) 4700 CONTINUE C TRANS = -2.0*TRANS C C WRITE TRANSITION PROBABILITIES: WRITE (6,13000) WRITE (6,37000) (JOTA(JAY),JAY=1,K) DO 4800 I1 = 1,K DO 4800 I2 = 1,K WRITE (6,38000) I1, I2, (NT(I1,I2,J),J=1,K) 4800 CONTINUE WRITE (6,39000) WRITE (6,37000) (JOTA(JAY),JAY=1,K) DO 4900 I1=1,K DO 4900 I2=1,K WRITE (6,40000) I1, I2, (TP(I1,I2,J),J=1,K) 4900 CONTINUE WRITE (6,63000) (PROB(IG),IG=1,K) WRITE (6,32000) TRANS WRITE (6,50000) WRITE (6,22000) XMN2LL C C ACCOUNT FOR LABELS OF BORDER OBSERVATIONS: FIRST = 0.0 DO 5000 J = 1,NC LABEL1 = ICLUS(1,J) PROBAB = PROB(LABEL1) FIRST = FIRST + ALOG(PROBAB) 5000 CONTINUE DO 5100 I = 2,NR LABEL1 = ICLUS(I,1) PROBAB = PROB(LABEL1) FIRST = FIRST + ALOG(PROBAB) 5100 CONTINUE FIRST = -2.0*FIRST C C C COMPUTE MODEL SELECTION CRITERIA: C COMPUTE VALUES CORRESPONDING TO NEW SEGMENTATION AND OLD C COVARIANCE MATRIX (OMIT ON FOR FIRST ITERATION) IF (ITER .EQ. 1) GO TO 5200 C FOR MODEL WITH COMMON COVARIANCE MATRIX C (NOT OPTIMIZED IN THIS PROGRAM; HOWEVER, IT IS CLEAR C THAT ONE SHOULD USE THE COMMON-COVARIANCE-MATRIX MODEL IF C THE VALUES OF THE MODEL SELECTION CRITERIA HERE FOR THAT MODEL C ARE LESS THAN THOSE FOR THE MODEL WITH DIFFERENT COVARIANCE C MATRICES): C AICOLD = XMN2OL + TRANS + FIRST + 2.0*NOPARM SCHOLD = XMN2OL + TRANS + FIRST + ALOG(XN)*NOPARM WRITE (6,59000) NOPARM WRITE (6,46000) AICOLD WRITE (6,47000) SCHOLD 5200 CONTINUE AIC = XMN2LL + TRANS + FIRST + 2.0*NOPARM WRITE (6,48000) AIC SCH = XMN2LL + TRANS + FIRST + ALOG(XN)*NOPARM WRITE (6,49000) SCH C C C COMPUTE PRECISION MATRICES FOR ALL GROUPS. C TERM=0.0 DO 5700 L=1,K IF (NG(L) .GT. IP) GO TO 5300 WRITE (6,52000) L,NG(L) GO TO 5600 5300 CONTINUE C DO 5400 IVAR=1,IP DO 5400 JV = 1,IP A(IVAR,JV) = SSD(L,IVAR,JV)/NG(L) 5400 CONTINUE IDET = 1 NRS1 = 0 CALL MATEQ(A,IP,20,JFLG,DET,IDET,IV,NRS1,P,20) DO 5500 IVAR=1,IP DO 5500 JV=1,IP PG(L,IVAR,JV) = P(IVAR,JV) 5500 CONTINUE C ET(L) = DLOG(DET*10**IDET) C 5600 CONTINUE C TERM = TERM + NG(L)*ET(L) 5700 CONTINUE C COMPUTE MODEL SELECTION CRITERIA WITH DIFFERENT COVARIANCE C MATRICES: C C PARAMETERS: C K MEAN VECTORS OF DIMENSION P AND K P-BY-P COVARIANCE MATRICES, C WHERE P IS THE NUMBER OF CHANNELS XM2LLD = N*IP*ALOG(2*PI) + N*IP + TERM XM2LLD = XM2LLD + TRANS + FIRST SCHD = XM2LLD + ALOG(XN)*NPRMDF AICD = XM2LLD + 2.0*NPRMDF WRITE (6,23000) WRITE (6,59000) NPRMDF WRITE (6,24000) XM2LLD WRITE (6,57000) AICD WRITE (6,58000) SCHD C C WRITE (6,25000) ITER DO 5800 IG = 1,K WRITE (6,21000) IG, (XMEAN(IG,IVAR),IVAR=1,IP) 5800 CONTINUE ITER = ITER + 1 IF (ITER .GE. 21) GO TO 5900 C UNLESS 20 ITERATIONS, HAVE ALREADY BEEN PERFORMED, GO BACK C AND DO ANOTHER. GO TO 900 5900 WRITE (6,44000) GO TO 6500 C C 6000 CONTINUE C C OUTPUT TO BE WRITTEN UPON CONVERGENCE: C WRITE (6,41000) ITER WRITE (6,30000) C COMPUTE AND WRITE COMMON COVARIANCE MATRIX C DO 6100 IVAR = 1,IP DO 6100 JV = 1,IP WGMS(IVAR,JV) = WGSS(IVAR,JV)/(N-K) 6100 CONTINUE DO 6200 IVAR=1,IP WRITE (6,31000) (WGMS(IVAR,JV), JV=1,IP) 6200 CONTINUE C WRITE (6,62000) C COMPUTE AND WRITE CLASS COVARIANCE MATRICES C DO 6400 L = 1,K DO 6300 IVAR=1,IP DO 6300 JV =1,IP WGMS(IVAR,JV) = SSD(L,IVAR,JV)/(NG(L)-1) 6300 CONTINUE WRITE (6,61000) L DO 6400 IVAR = 1,IP WRITE (6,31000) (WGMS(IVAR,JV),JV=1,IP) 6400 CONTINUE C C C C C C WRITE (6,45000) K C 6500 CONTINUE C WRITE (6,42000) 6600 STOP C C 10000 FORMAT(//1X,'SEGMENTATION:'/) 11000 FORMAT(/,1X,'ROW: COLUMN: '/, 4X, (40I3/) ) 12000 FORMAT(1X, I3, (40I3/) ) 13000 FORMAT(//1X,'TRANSITIONS') 14000 FORMAT(/1X,'NUMBER OF CHANNELS = ',I2/) 15000 FORMAT(2X,I2) 16000 FORMAT(3X,I3) 17000 FORMAT('1','###################################################', X//,1X,'PROGRAM IMSEG.DET.AUTO '/ X,1X,'FOR IMAGE SEGMENTATION '/ X,1X,'USING DISTANCE IN THE METRICS OF THE COVARIANCE MATRICES'/ X,1X,'ADJUSTED BY THE DETERMINANTS '/ X//,1X,'DEVELOPED AND PROGRAMMED BY '// X1X,' DR. STANLEY L. SCLOVE 312/996-2681'/ X1X,' DEPARTMENT OF QUANTITATIVE METHODS 312/996-2676'/ X1X,' UNIVERSITY OF ILLINOIS AT CHICAGO '/ X1X,' BOX 4348, CHICAGO, IL 60680 '// X//,1X,'VERSION 1.6 17-SEP-82 '//) 18000 FORMAT('1',//,1X,'MIN AND MAX NUMBER OF CLASSES TO BE TRIED ARE ', XI2,' AND ',I2//) 19000 FORMAT(1H1,'K = ',I2,' CLASSES'//) 20000 FORMAT(/1X,'INITIAL MEAN VECTOR FOR CLASS ',I2,': ',(4E14.5/)) 21000 FORMAT(1X,'MEAN VECTOR FOR CLASS ',I2,': ',(8E13.5/)) 22000 FORMAT(/,1X,'MINUS 2 LOG LIKELIHOOD FOR MODEL WITH COMMON', X' COVARIANCE MATRIX= ', E13.5//) 23000 FORMAT(//1X,'FOR MODEL WITH DIFFERENT COVARIANCE MATRICES:'//) 24000 FORMAT(/,1X,'MINUS 2 LOG LIKELIHOOD FOR MODEL WITH DIFFERENT ', X'COVARIANCE MATRICES = ', E13.5//) 25000 FORMAT(///,1X,'ITERATION ', I2,//) 26000 FORMAT(1X,'MEAN VECTOR FOR CLASS ',I2,': ',(8E13.5/)) 27000 FORMAT(/,1X,'NUMBERS IN CLASSES:'/,(9I12/)/) 28000 FORMAT(18A4) 29000 FORMAT(3X,I2) 30000 FORMAT(///,1X,'COMMON COVARIANCE MATRIX (DIVISOR IS DF):',//) 31000 FORMAT(1X,(8E13.5/)) 32000 FORMAT(//,' CONTRIBUTION OF TRANS. PROBS. TO LOG LIKELIHOOD =', XE15.5/) 33000 FORMAT(1X,'NUMBER OF ROWS = ',I4) 34000 FORMAT(1X,'NUMBER OF COLUMNS = ',I3/) 35000 FORMAT(1X,18A4) 36000 FORMAT(/9X,(9I7/)) 37000 FORMAT(/9X,(9I7/)) 38000 FORMAT(1X, 2I4, (9I7/)) 39000 FORMAT(//1X,'TRANSITION PROBABILITIES') 40000 FORMAT(1X,2I4,3X,(9F7.4/)) 41000 FORMAT(/1X,'CONVERGENCE: NO CASE CHANGED CLASSES AFTER ', X'ITERATION ',I2,'. SOME ADDITIONAL RESULTS ARE PRINTED BELOW.'//) 42000 FORMAT(/1X,'PROGRAM ENDED SUCCESSFULLY.') 43000 FORMAT(/,1X,'JFLG = ',I2,'. IF JFLG=0, COMPUTATION OF DET', X' WENT WELL; OTHERWISE, THERE WAS TROUBLE OR MATRIX WAS ', X'ILL-CONDITIONED.'//) 44000 FORMAT(1X,'ROUTINE HAS NOT CONVERGED IN 20 ITERATIONS. STOP') 45000 FORMAT(/,1X,'PROGRAM ENDED SUCCESSFULLY FOR THE CASE ', X'K = ',I2,'.'/) 46000 FORMAT(/1X,'AIC WITH NEW LABELS AND OLD DISTRIBUTIONAL ', X'PARAMETERS: ', E15.5) 47000 FORMAT(/1X,'SCHWARZ CRITERION WITH NEW LABELS AND OLD ', X'DISTRIBUTIONAL PARAMETERS: ', E15.5) 48000 FORMAT(/1X,'AIC WITH NEW LABELS AND NEW DISTRIBUTIONAL ', X'PARAMETERS: ',E15.5) 49000 FORMAT(/1X,'SCHWARZ CRITERION WITH NEW LABELS AND NEW ', X'DISTRIBUTIONAL PARAMETERS: ',E15.5) 50000 FORMAT(//1X,'FOR MODEL WITH COMMON COVARIANCE MATRIX ', X'(NOT OPTIMIZED IN THIS PROGRAM): '//) 51000 FORMAT(1X,'AIC FOR MODEL WITH COMMON COVARIANCE MATRIX = ',E15.5/) 52000 FORMAT(/1X,'CLASS ',I2,' CONTAINS ONLY ',I4,' OBSERVATIONS: ', X'PRECISION MATRIX FROM PREVIOUS ITERATION IS BEING RETAINED'/) 53000 FORMAT(/1X,'MODEL WITH COMMON COVARIANCE MATRIX IS NOT'/ X1X,'OPTIMIZED IN THIS PROGRAM; HOWEVER, IT IS CLEAR'/ X1X,'THAT ONE SHOULD USE IT IF HERE THE VALUES OF THE MODEL-'/ X1X,'SELECTION CRITERIA FOR THAT MODEL ARE LESS THAN THOSE '/ X1X,'FOR THE MODEL WITH DIFFERENT COVARIANCE MATRICES,'/ X1X,'WHICH IS OPTIMIZED HERE.'//) 54000 FORMAT(/1X,'MINIMUM FOR EACH CHANNEL: ',/) 55000 FORMAT(/1X,'MAXIMUM FOR EACH CHANNEL: ',/) 56000 FORMAT(1X,'SCHWARZ CRITERION FOR MODEL WITH COMMON COVARIANCE ', X'MATRIX = ',E15.5/) 57000 FORMAT(1X,'AIC FOR MODEL WITH DIFFERENT COVARIANCE MATRICES = ', XE15.5/) 58000 FORMAT(1X,'SCHWARZ CRITERION FOR MODEL WITH DIFFERENT ', X'COVARIANCE MATRICES = ', E15.5/) 59000 FORMAT(/,1X,'NUMBER OF PARAMETERS = ',I4//) 60000 FORMAT(1X,'NO OBSERVATIONS IN GROUP ',I3,'. STOP') 61000 FORMAT(//1X,'COVARIANCE MATRIX FOR CLASS ',I2/) 62000 FORMAT(///1X,'CLASS COVARIANCE MATRICES ', X' (DIVISORS ARE ONE LESS THAN NUMBER IN GROUP):'/) 63000 FORMAT(/1X,'MARGINAL PROB.VECTOR:',(9F11.4/)) END SUBROUTINE MATEQ(A,M,N,JFLG,DET,IDET,IV,NRS1,P,LL) C SUBROUTINE MATEQ IS DMATEQ FROM THE UICC SUBROUTINE LIBRARY. C C SUBROUTINE DMATEQ C ***************** C THIS ROUTINE WILL SOLVE A REAL*8 SYSTEM OF LINEAR EQUATIONS,COMPUTE C THE DETERMINANT, WITHOUT UNDERFLOW OR OVERFLOW, OF A REAL*8 MATRIX, C AND/OR INVERT A REAL*8 MATRIX. C CALLING SEQUENCE: C CALL DMATEQ(A,N,IA,JFLG,DET,IDET,IV,NRS,P,IP) WHERE; C A (INPUT) - IS THE REAL*8 MATRIX ON WHICH THE ROUTINE IS C TO WORK. IN THE PROCESS OF COMPUTATION THE C CONTENTS OF THIS MATRIX ARE DESTROYED. C N (INPUT) - IS AN INTEGER*4 VARIABLE WHICH SPECIFIES THE C ORDER OF THE A MATRIX. C IA (INPUT) - IS AN INTEGER*4 VARIABLE WHICH SPECIFIES THE C ACTUAL ROW DIMENSION OF A AS DIMENSIONED IN C THE CALLING PROGRAM. IA MUST BE GREATER THAN C OR EQUAL TO N. C JFLG (OUTPUT) - IS AN INTEGER*4 RETURN CODE VARIABLE. UPON C RETURN FROM DMATEQ IF; C JFLG=0, ALL WENT WELL. C JFLG=1, THE A MATRIX WAS SINGULAR OR NEAR C SINGULAR AND THE COMPUTATIONS COULD NOT BE C COMPLETED. THE CONTENTS OF THE VARIABLES C A, DET, IDET AND P ARE MEANINGLESS. C DET (OUTPUT) - IS A REAL*8 VARIABLE WHICH CONTAINS THE C DETERMINANT OF A. (SEE IDET) C IDET (INPUT) - IS AN INTEGER*4 VARIABLE. ON INPUT IF; C IDET=0, NO DETERMINANT IS CALCULATED. C IDET NOT 0, THE DETERMINANT OF A IS COMPUTED. C ON OUTPUT IDET CONTAINS THE POWER OF 10 C THAT DET SHOULD BE MULTIPLIED BY TO GIVE THE C CORRECT VALUE OF THE DETERMINANT. I.E. C DET(A)=DET*10.0D0**IDET. C IF DET(A) CAN BE COMPUTED WITHOUT UNDER OR C OVERFLOW, THEN IDET=0 OTHERWISE IDET IS SET C TO THE PROPER VALUE SO THAT NO UNDER OR OVER- C FLOW WILL OCCUR IN COMPUTING DET. C IV (INPUT) - IS AN INTEGER*4 WORK ARRAY WHICH SHOULD BE C DIMENSIONED AT LEAST IV(N). C C NRS (INPUT) - IS AN INTEGER*4 VARIABLE WITH THE FOLLOWING C INTERPRETATION: C NRS>0, SOLVE A SYSTEM OF LINEAR EQUATIONS C WITH NRS RIGHT HAND SIDES. C NRS=0, INVERT THE A MATRIX. C NRS<0, ONLY COMPUTE THE DETERMINANT OF A. C IN THIS CASE IDET MUST BE DIFFERENT FROM 0. C P (INPUT) - IS A REAL*8 ARRAY WITH THE FOLLOWING INTER- C PRETATION: C IF NRS>0, THEN P CONTAINS THE NRS RIGHT HAND C SIDES STORED BY COLUMNS. IN THIS CASE P MUST C BE DIMENSIONED AT LEAST P(N,NRS). ON RETURN C THE COLUMNS OF P ARE REPLACED BY THE RESPEC- C TIVE SOLUTIONS. C IF NRS=0, THEN P MUST BE DIMENSIONED AT LEAST C P(N,N). ON RETURN P WILL CONTAIN THE INVERSE C OF A. C IF NRS<0,THEN P NEED ONLY BE A DUMMY VARIABLE C IN THIS CASE P IS NEVER ACCESSED BY DMATEQ. C IP (INPUT) - IS AN INTEGER*4 VARIABLE WHICH CONTAINS THE C ACTUAL ROW DIMENSION OF P AS DIMENSIONED IN C THE CALLING PROGRAM. IP MUST BE GREATER THAN C OR EQUAL TO N. C NOTE: IMMEDIATELY ON RETURN FROM DMATEQ THE CONDITION CODE FLAG, C JFLG, SHOULD BE INTERROGATED. IF JFLG=1, THEN THE ROUTINE C COULD NOT COMPUTE A SOLUTION. C METHOD - THE ALGORITHM USED IS GAUSSIAN ELIMINATION WITH PARTIAL C -1 C PIVOTING. IN ESSENCE THE ROUTINE GENERATES A MATRIX L SUCH C -1 C THAT L *A = U, WHERE U IS AN UPPER TRIANGULAR MATRIX. THEN IT C SOLVES THE SYSTEM A*X = P BY MEANS OF THE EQUIVALENT SYSTEM C -1 -1 C U*X = L *A*X = L *P BY BACK SUBSTITUTION. C -1 C THE L MATRIX CAN BE WRITTEN AS A PRODUCT OF THE FORM C -1 C L = L *P *....*L *P WHERE EACH P IS A PERMUTATION C N-1 N-1 1 1 K C MATRIX OBTAINED BY INTERCHANGING AT MOST TWO ROWS OF THE C IDENTITY MATRIX. ( THIS REPRESENTS THE INTERCHANGING OF TWO C ROWS). THE L MATRICES ARE ELIMINATION MATRICES WHICH ARE C K C CHOSEN TO INTRODUCE ZEROS IN THE LAST N-K ENTRIES OF THE K-TH C COLUMN OF THE MATRIX. C C -1 -1 C THE CALCULATIONS OF L *A AND L *P ARE DONE BY PERFORMING C THE PERMUTATIONS ON A AND P RESPECTIVELY. THE ACTUAL L AND P C K K C ARE NOT COMPUTED. C SUBROUTINES CALLED: DMATDT C REFERENCE: C G. W. STEWART, INTRODUCTION TO MATRIX COMPUTATIONS, C ACADEMIC PRESS, 1973. REAL*8 A(N,1),DET,P(LL,1) REAL*8 DNORM,DEN,DMULT,DSUM,DISIGN DIMENSION IV(1) NRS=NRS1 IF (NRS.EQ.0) IDET=1 DISIGN=1.0D+00 DET=0.0D+00 JFLG=0 C C JFLG IS A TROUBLE FLAG.UPON EXIT IF JFLG=0 THEN THE MATRIX WAS PROCESS C WITHOUT TROUBLE.IF JFLG=1 EITHER THE MATRIX IS SINGULAR OR TROUBLE C OCCURED.ISIGN=-ISIGN EVERY TIME A ROW IS INTERCHANGED.THIS IS USED TO C INSURE THAT THE DETERMINANT HAS THE PROPER SIGN. C M1=M-1 DO 100 I=1,M 100 IV(I)=I IF (NRS) 500,200,500 200 DO 300 I=1,M DO 300 J=1,M 300 P(I,J)=0.0D+00 DO 400 I=1,M 400 P(I,I)=1.0D+00 NRS=M C C INSTEAD OF ACTUALLY INTERCHANGING ROWS A POINTER ARRAY IS USED TO KEEP C TRACK OF THE ROW POSITIONS. C C BEGIN ELIMINATION LOOP. C 500 DO 1200 K=1,M1 ICOL=K IPCOL=K C C SEARCHING FOR LARGEST ELEMENT IN ABSOLUTE VALUE IN COLUMN K. C DNORM=A(IV(K),K) IFLG=0 KK=K+1 DO 600 J=KK,M IF (DABS(A(IV(J),K)).LE.DABS(DNORM)) GO TO 600 IFLG=1 IPCOL=IV(J) DNORM=A(IPCOL,K) 600 CONTINUE C C IF IFLG=0 NO ROW INTERCHANGE TOOK PLACE.IF IFLG=1 A ROW INTERCHANGE C TOOK PLACE AND THE POINTER ARRAY IV MUST BE UPDATED. C IF (IFLG.EQ.0) GO TO 800 ISAVE=IV(ICOL) IV(ICOL)=IPCOL ICOL1=ICOL+1 DO 700 L=ICOL1,M IF (IV(L).EQ.IPCOL) IV(L)=ISAVE 700 CONTINUE DISIGN=-DISIGN 800 IF (DNORM.EQ.0.0D+00) GO TO 1900 C C BEGIN ELIMINATION OF ROW BELOW IV(K).DEN IS THE PIVOT ELEMENT. C K1=K+1 DO 1100 IM=K1,M C C BEFORE ACTUALLY ELIMINTING WE CHECK TO SEE IF A(IV(IM),K) HAS ALREADY C BEEN ANIHALATED. C IF (A(IV(IM),K).EQ.0.0D+00) GO TO 1100 C C CACULATE ELIMINATION FACTOR. C DMULT=-A(IV(IM),K) C C WE NOW CALCULATE VALUE OF OTHER ELEMENTS IN ROW IV(IM). C DO 900 NN=K1,M 900 A(IV(IM),NN)=(DMULT*A(IV(K),NN))/DNORM+A(IV(IM),NN) IF (NRS.LE.0) GO TO 1100 DO 1000 IN=1,NRS 1000 P(IV(IM),IN)=(DMULT*P(IV(K),IN))/DNORM+P(IV(IM),IN) 1100 CONTINUE 1200 CONTINUE C C CALCULATE VALUE OF DETERMINANT. C IF (A(IV(M),M).EQ.0.0D0) GO TO 1900 DET=DISIGN IF (IDET.NE.0) CALL DMATDT(A,N,M,DET,IV,IDET) IF (DET.EQ.0.0D+00) GO TO 1900 IF (NRS.LE.0) GO TO 2000 C C WE START SOLVING RIGHT HAND SIDES.THE SOLUTION REPLACES THE RIGHT HAND C VECTOR. C 1300 N1=M-1 DO 1600 JJ=1,NRS C C BEGIN BACK SUBSTITUTION. C P(IV(M),JJ)=P(IV(M),JJ)/A(IV(M),M) DO 1500 I=1,N1 DSUM=0.0D+00 DO 1400 J=1,I 1400 DSUM=DSUM-A(IV(M-I),M-J+1)*P(IV(M-J+1),JJ) 1500 P(IV(M-I),JJ)=(P(IV(M-I),JJ)+DSUM)/A(IV(M-I),M-I) 1600 CONTINUE DO 1800 JJ=1,NRS DO 1700 IND=1,M 1700 A(IND,1)=P(IV(IND),JJ) DO 1800 IND=1,M 1800 P(IND,JJ)=A(IND,1) RETURN 1900 JFLG=1 IDET=0 2000 RETURN END SUBROUTINE MATDT(A,IA,N,DET,IV,IDET) C SUBROUTINE MATDT IS DMATDT FROM THE UICC SUBROUTINE LIBRARY. REAL*8 A(IA,1),DET,B,LOG16 INTEGER*4 IV(1),K EQUIVALENCE (B,K) NUM=16777216 LOG16=.120411998265592457D+01 IF (A(IV(N),N).EQ.0.0D+00) GO TO 300 L=0 DO 100 I=1,N B=DABS(A(IV(I),I)) K=K/NUM-64 L=L+K 100 DET=DET*(A(IV(I),I)/16.0D+00**K) B=DABS(DET) K=K/NUM-64 IW=L+K IF ((IW.LT.-64).OR.(IW.GT.63)) GO TO 200 DET=DET*16.0D+00**L IDET=0 GO TO 400 200 DET=DET*16.0D+00**(-K) IDET=L+K B=IDET*LOG16 IDET=B B=B-DFLOAT(IDET) DET=DET*1.0D+01**B GO TO 400 300 DET=0.0D+00 IDET=0 400 RETURN END