C	PROGRAM DYN1.FT
C	----------------
C
C	###SUBROUTINE   DYN1
C
C
C
C	PETER LEMKIN
C	IMAGE PROCESSING UNIT, DCBD
C	NATIONAL CANCER INSTITUTE
C	NATIONAL INSTITUTES OF HEALTH
C	9000 ROCKVILLE PIKE
C	BETHESDA, MD. 20014
C
C	MAY 5, 1978
C	MAY 4, 1978
C	MAY 2, 1978
C	MAY 1, 1978
C	APRIL 28, 1978
C	APRIL 27, 1978
C	APRIL 20, 1978
C	APRIL 18, 1978
C	APRIL 13, 1978
C
C
C
C	PURPOSE
C	-------
C	COMPUTE THE OPTIMAL THRESHOLD RETURNED IN "ISDEV"
C	FOR THE WINDOW AT
C	(KX1:KX2,KY1:KY2) FOR (IBM1,IBYTE,IX,IY).
C
C	PARAMETERS PASSED TO DYN1 FROM DYNBND
C	-------------------------------------
C	MODEPTILE 	- %TILE THRESHOLD
C	MODEDPTILE 	- DELTA %TILE ALLOWABLE ERROR.
C	MODECOUNT	- POSITIVE COUNT USED TO RESET LSPALCOUNTER.
C	LSPALCOUNTER	- NEGATIVE VALUED UPCOUNTER TO RECOMPUTE
C			  HISTOGRAM WHEN REACH ZERO.
C
C	OPDEFS
C	------
S	OPDEF	HPL	6360
S	OPDEF	VPL	6362
S	OPDEF	HSL	6361
S	OPDEF	VSL	6363
C
S	OPDEF	LDXP	6443
S	OPDEF	LDYP	6444
C
S	OPDEF	DISP1	6435
S	OPDEF	DISP2	6436
C
S	OPDEF	TADI	1400
S	OPDEF	ISZI	2400
S	OPDEF	DCAI	3400
C
C
S	OPDEF	SWAB	7431
S	OPDEF	SWBA	7447
S	OPDEF	DAD	7443
S	OPDEF	DST	7445
S	OPDEF	CLAMQ	7621
S	OPDEF	MUY	7405
S	OPDEF	DVI	7407
C
S	OPDEF	LSR	7417
S	OPDEF	MQA	7501
S	OPDEF	MQL	7421
S	OPDEF	BSW	7002
S	OPDEF	DVI	7407
C
C
C
C	[0] ENTRY
S	ENTRY DYN1
S	CPAGE 2
S DYN1,	BLOCK 2
C
C
C	[1] DO HIST OF WINDOW ==>BM3
C	COPY /H SWITCH
	ISWH=ISW(8)
	JXSAVE=IX
	JYSAVE=IY
C
C
C	[1.1] TEST IF RECOMPUTE THE HISTOGRAM AND ITS PARAMETERS;
C	IF (LSPALCOUNTER_LSPALCOUNTER+1)=0
C		THEN RESET COUNTER AND COMPUTE HISTOGRAM;
S	ISZ \LSPALCOUNTER
S	JMP \112 /CONTINUE
C
C
C	[1.1.1] COMPUTE HISTOGRAM AND ITS PARAMETERS
S	JMS DOHIST
C
C	LSPALCOUNTER_(IF KDEVOUT THEN -MODECOUNT ELSE -1);
S	TAD \KDEVOUT
S	SZA CLA
S	TAD \MODECOUNT
S	SNA
S	IAC /1
S	CIA
S	DCA \LSPALCOUNTER
C
C
C	[1.2] IF /H THEN DRAW THE HISTOGRAM
S	TAD \ISWH
S	SNA CLA
S	JMP \112 /GOTO [1.1.2]
C
C	FIRST ZERO BM3
	CALL BMOMNI(3,0, 0,0,0, IBUF1,4)
C	COMPUTE:	CURSYM=252.0/FLOAT(IH(MODE+1))
S	TAD \MODE
S	TAD PIH
S	DCA 7
S	CPAGE 2
S	6211
S	TADI 7
S	DCA \IZ
	CURSYM=252.0/FLOAT(IZ)
C
	DO 431 IX1=1,256
	IX=IX1-1
C	COMPUTE:	IY2=IH(IX1)
S	TAD PIH
S	CPAGE 5
S	TAD \IX
S	DCA 7
S	TADI 7
S	DCA \IY2
C
C	IF IY2=0
C		THEN NOP;
S	TAD \IY2
S	SNA CLA
S	JMP \431 /NOP
C
	IY2=CURSYM*FLOAT(IY2+4)
	IZ=255
	MEM=3
	IBYTE=0
C
	DO 431 IY1=4,IY2
	IY=257-IY1
	CALL PACK2D
431	CONTINUE
C
C
C	[1.2] INIT PARAMETERS
C	COMPUTE THE RELATIVE WINDOW SIZE
112	FC=FLOAT(NTOT)/100.0
	MODEMEDIAN=(MING+MAXG)/2
C
C	ESTIMATE SLICE THRESHOLD AS MEDIAN
	ITSLICE=MODEMEDIAN
C
C
C	[1.3] CONTINUE USING LAST VALUE OF ESTIMATED ITSLICE
C	AND HISTOGRAM PARAMETERS
C	CLEAR THE "CONSTRAINT FAILED" VARIABLE
	IOPTOP=0
C
C	DEFINE MAX NEG. COUNT FOR HEURISTIC SEARCH FOR OPT. THRESHOLD.
	MAXLOOP=-(MAXG-MING)/2
C
C
C
C	[2] TEST IF CAUGHT IN LOOP OF CONSTRAINT FAILURE.
349	CONTINUE
C	IF (MAXLOOP_MAXLOOP+1)=0
C		THEN EXIT
S	ISZ \MAXLOOP
S	JMP \1349 /CONTINUE, NOT YET
C
C	FATAL EXIT - NOTE: THRESHOLD IS NEAR OPTIMUM
	IOPTOP=8
	GOTO 360
C
S	CPAGE 3
S\1349,	JMS TTYCTL
S	JMP \998
C
C
C	[2.1] COMPUTE CURRENT PTILE AT ITSLICE
	IVAL=0
	DO 339 IP=MING,ITSLICE
S	TAD PIH
S	CPAGE 5
S	TAD \IP
S	DCA 7
S	TADI 7
S	TAD \IVAL
S	DCA \IVAL
339	CONTINUE
C	COMPUTE THE PTILE
	IP=FLOAT(IVAL)/FC
C
C
C	[2.2] COMPUTE SLICED NGH I10[0:8] AS BINARY NGH.
C	RESTORE THE STATE
	IX=JXSAVE
	IY=JYSAVE
	MEM=IBM1
	IBYTE=IHGH1
C	COMPUTE BINARY NGH FROM (BMI SLICE ITSLICE)
	CALL GETI1
	IA=0
S	TAD \ITSLICE
S	CIA
S	DCA 25 /SAVE IT
C
	DO 348 IZ=1,9
S	CLA CMA
S	TAD PI10
S	CPAGE 12
S	TAD \IZ
S	DCA 7
S	TAD 25 /-ITSLICE
S	TADI 7 /I10(IZ)
S	SMA CLA
S	IAC /> THR
S	DCAI 7 /I10(IZ)
S	TADI 7 /I10(IZ)
S	TAD \IA
S	DCA \IA
348	CONTINUE
C
C
C	[2.3] TEST CONSTRAINTS...
C
C	[2.3.1] TEST #1
C	IF (IP < (MODEPTILE-MODEDPTILE))
C		THEN ITSLICE_ITSLICE+1 GOTO [2]
S	TAD \MODEDPTILE
S	CIA
S	TAD \MODEPTILE
S	CIA
S	TAD \IP
S	SMA CLA
S	JMP \345 /NO
	ITSLICE=ITSLICE+1
	IOPTOP=1
	GOTO 358
C
C
C	[2.3.2] TEST #2
C	IF (IP > MODEPTILE+MODEDPTILE)
C		THEN ITSLICE_ITSLICE-1 GOTO [2]
S\345,	TAD \MODEPTILE
S	TAD \MODEDPTILE
S	CIA
S	TAD \IP
S	SPA CLA
S	JMP \346 /NO
	ITSLICE=ITSLICE-1
	IOPTOP=2
	GOTO 358
C
C
C	[2.3.3] TEST #3
C	IF (I10+I11+I12+I13+I14+I15+I17)=0
C		THEN ITSLICE_ITSLICE-1;
C		     GOTO [2];
346	IZ=I10+I11+I12+I13+I14+I15+I16+I17
S	TAD \IZ
S	SZA CLA
S	JMP \247 /OK
	ITSLICE=ITSLICE-1
	IOPTOP=3
	GOTO 358
C
C
C	[2.3.4] TEST #4
C	IF IA=9
C		THEN ITSLICE_ITSLICE+1;
C		     GOTO [2];
S\247,	TAD \IA
S	TAD (-D9
S	SZA CLA
S	JMP \248 /OK
	ITSLICE=ITSLICE+1
	IOPTOP=4
	GOTO 358
C
C
C	[2.3.5] TEST #5
C	IF I18=0
C		THEN ITSLICE_ITSLICE-1, GOTO [2];
S \248,	TAD \I18
S	SZA CLA
S	JMP \359 /OK
	ITSLICE=ITSLICE-1
	IOPTOP=5
	GOTO 358
C
C
C	[2.3.8] FAILED! GOTO [2] TO TRY AGAIN.

358	CONTINUE
C****DEBUG*****
S	6344 /FBW4
S	AND (4000
S	SNA CLA
S	JMP \349 /NO
C	PRINT ON TTY:. INCREMENTAL INFO
	LCNT=1
S	JMS PRINTHSWITCH
C*****************
	GOTO 349
C
C
C	[2.3.9] SUCCESSFUL!
359	CONTINUE
	IOPTOP=9
C	NOTE: I10:I18 HAS BINARY IMAGE.
C
C
C	[2.2] IF /H THEN PRINT THE RESULTS
360	LCNT=3
S	TAD \ISWH
S	SZA CLA
S	JMS PRINTHSWITCH
C
C
C	[3] IF /H THEN DO SLICE OF WINDOW ==>BM3
S	TAD \ISWH
S	SNA CLA
S	JMP \998 /NO SLICE
C
	DO 320 IY1=KY1,KY2
	IY=IY1-1
	DO 320 IX1=KX1,KX2
	IX=IX1-1
	MEM=IBM1
	IBYTE=IHGH1
	CALL FETCH2D
C
C	IF IZ < ITSLICE
C		THEN IZ<=0;
S	TAD \ITSLICE
S	CIA
S	TAD \IZ
S	SMA CLA
S	TAD (D255
S	DCA \IZ
	MEM=3
	IBYTE=0
	CALL PACK2D
	MEM=IBM1
	IBYTE=IHGH1
320	CONTINUE
C
C
C	[4] CONTINUE
998	IX=JXSAVE
	IY=JYSAVE
C	RETURN THE THRESHOLD
	ISDEV=ITSLICE
S	RETRN DYN1
C	************************************************************
C	SUBROUTINE:	T T Y C T L	(INTERNAL)
C	************************************************************
C
S	CPAGE 3
S RTTYC,	JMP I TTYCTL
S TTYCTL,	0000	/ENTRY
C
S	KSF	/ANYTHING TYPED
S	JMP NORMAL	/NO, RETURN NORMALLY
S	KRB	/GET TYPED CHARACTER
S	AND (0177	/TAKE CARE OF PARITY PROBLEMS
S	TAD (-17	/TEST FOR CTRL/O
S	SNA	/SKIP IF NOT CTRL/O
S	JMP RTTYC	/ABORT CALLING ROUTINE (ERROR RETURN)
S	TAD (-4	/TEST FOR CTRL/S [-17-4=-23(OCTAL)]
S	SZA CLA	/SKIP IF CTRL/S
S	JMP NORMAL	/NOT CTRL/O OR CTRL/S SO RETURN NORMALLY
C
S SLEEP,KSF	/WAIT FOR CTRL/Q
S	JMP SLEEP	/KEEP WAITING
S	KRB	/READ CHARACTER
S	AND (0177
S	TAD (-17	/IS IT A CTRL/O?
S	SNA	/SKIP IF NOT
S	JMP RTTYC	/YES, ABORT
S	TAD (-2	/TEST FOR CTRL/Q (-17-2=-21 OCTAL)
S	SZA CLA	/SKIP IF SO
S	JMP SLEEP	/NOPE, KEEP SLEEPING
C
S NORMAL,INC TTYCTL	/INCREMENT RETURN ADDRESS FOR NORMAL RETURN
S	CLA	/SAFETY VALVE
S	JMP RTTYC	/RETURN
C	******************************************************
C	*SUBROUTINE     D O H I S T
C	*****************************************************
C
C	PURPOSE
C	-------
C	COMPUTE
C	VARIABLES
C	---------
C	IH[0:255]	-SMOOTHED HISTOGRAM OF PIXEL DATA IN WINDOW.
C	MING		- LUB(IH)
C	MAXG		- GLB(IH)
C
C	MODE		- LARGEST PEAK OF IH.
C	MV		- IH[MODE]
C	STDDB		- 1-SIDED STD DEV FROM LEFT OF MODE.
C
C	MEANO		- MEAN OF ALL IH
C	STDDO		- 2-SIDED STD DEV OF ALL IH FROM MEANO.
C
C	NTOT		- TOTAL # SMOOTHED PIXELS
C
C
C
C
C	OF THE GRAY VALUES IN IBM1 DEFINED BY (KX1:KX2,KY1:KY2)
C	HISTOGRAM. RETURN THE VALUES IN
C
C	[DH.0] ENTRY
S	CPAGE 3
S RDOHIST,	JMP I DOHIST
S DOHIST,	0
C
C
C	[DH.0.1] INIT
	MEM=IBM1
	IBYTE=IHGH1
	NTOT=0
	MV=0
	MING=2047
	MAXG=0
	MODE=0
C	ZERO SUM IZ*IH[IZ]
	IC=0
S	DCA \IC#
C
C	ZERO ==>IH[1:256]
	DO 100 IZ=1,256
S	CLA CMA
S	TAD PIH
S	CPAGE 5
S	TAD \IZ
S	DCA 7
S	DCAI 7
100	CONTINUE
C
C
C	[DH.1] LINE LOOP
	DO 200 IY1=KY1,KY2
	IY=IY1-1
S	TAD \IY
S	DISP2
C
C
C	[DH.2] PIXEL LOOP
	DO 200 IX1=KX1,KX2
	IX=IX1-1
S	TAD \IX
S	DISP1
	CALL FETCH2D
C
C	COMPUTE: IH(IZ)=IH(IZ)+1
S	TAD PIH
S	CPAGE 7
S	TAD \IZ
S	DCA 7
S	ISZI 7
S	NOP
S	TADI 7
S	DCA \IA /SAVE IH(IZ)
200	CONTINUE
C
C
C	[DH.2.1] SMOOTH H[.] AND THEN COMPUTE MING, MAXG, MODE.
	I11=0
	DO 210 IZ=2,253
C
C
C	[DH.2.1.1] IF /A
C		 THEN SMOOTH H[.] USING WEIGHTED AVERAGE USING
C			H'[I]=(H[I-1]+2H[I]+H[I+1])/4
C	FETCH H[I], H[I+1]
S	CLA CMA
S	TAD PIH
S	TAD \IZ
S	DCA 11
S	CPAGE 10
S	6211
S	TADI 11
S	DCA \I12
S	TAD 11
S	DCA 7
S	TADI 11
S	DCA \I13
C
C	IF /A
C		THEN AVG
C		ELSE IA_I12
	IA=I12
S	TAD \ISW
S	SNA CLA
S	JMP STOREIT /NO AVG
C
	IA=I11+(I12+I12)+I13
S	TAD \IA
S	CLL RAR ;CLL RAR /DIVIDE BY 4
S	DCA \IA
C
S	CPAGE 4
SSTOREIT,	TAD \IA
S	DCAI 7 /IH[IZ+1]
C
C	ROTATE LEFT VALUES
C	H[I-1]<==H[I].
	I11=I12
C
C	GET TOTAL NUMBER OF PIXELS
	NTOT=NTOT+IA
C
C
C	[DH.2.1.2] IF IA NEQ 0
C			 THEN COMPUTE MING, MAXG, MODE.
S	TAD \IA
S	SNA CLA
S	JMP \210 /NOP
C
C	COMPUTE: IC[1:2]=SUM IZ*IH[IZ] (I.E. IA*IA)
S	CLA CMA
S	CPAGE 12
S	TAD \IZ
S	SWAB
S	MUY
S	\IA
S	DAD
S	\IC
S	DST
S	\IC
S	CLAMQ
S	SWBA
C
C	COMPUTE: MING=MIN(MING, IZ)
C	         MAXG=MAX(MAXG,IZ).
S	TAD \IZ
S	CIA
S	TAD \MING
S	SPA CLA
S	JMP \202 /"MAX"
C
C	IS MIN
	MING=IZ-1
	GOTO 203
C
C	IS MAX
S\202,	TAD \MAXG
S	CIA
S	TAD \IZ
S	SPA CLA
S	JMP \203 /NEITHER
	MAXG=IZ-1
C
C	COMPUTE: MODE= (G|H(G) IS MAX OVER ENTIRE RANGE);
C	IF H(IZ) > MV
C		THEN MV_H(IZ), M1_IZ;
S\203,	TAD \IA /H(IZ)
S	CIA
S	TAD \MV
S	SMA CLA
S	JMP \210 /NO
	MV=IA
	MODE=IZ-1
210	CONTINUE
C
C
C	[DH.2.2] COMPUTE: MEANO=(1/NTOT)*SUM(IZ*IH[IZ])
C			       = (1/NTOT)*IC[1:2];
	CALL DPCVRT(IC,CURSYM,-1)
	MEANO=CURSYM/FLOAT(NTOT)
C	[DH.3] DONE COMPUTING IH. ESTIMATE 1-SIDED STD DEV FOR BACKGROUND
C	COMPUTE: VARIANCE=(1/N)*SUM(MV-IH(I))**2 FOR I=MING:MODE.
	IB=0
S	DCA \IB#
C
C	COMPUTE OVER 1:256 SO NO PROBLEM W/DO LOOP.
	I11=MING+1
	I12=MODE+1
	DO 310 I13=I11,I12
	IZ=I13-1
S	TAD PIH
S	CPAGE 5
S	TAD \IZ
S	DCA 7
S	TADI 7
S	CIA
S	TAD \MV
S	DCA \IA /CONTAINS:   (MV-IH(I13-1)).
C
C	COMPUTE: IA=IA**2=IA*IA
S	CPAGE 12
S	TAD \IA
S	SWAB /WITH MQL
S	MUY
S	\IA /CONTAINS (MV-IH(I13-1))
S	DAD  /NOW COMPUTE: IB[1:2]=IB[1:2] + (MV-IH(I13-1))**2
S	\IB
S	DST
S	\IB
S	CLAMQ
310	CONTINUE
S	SWBA
	N=MODE-MING+1
C
C
C	[DH.3.1] VARIANCE=IB[1:2]/N
	CALL DPCVRT(IB,CURSYM,-1)
	CURSYM=CURSYM/FLOAT(N)
C
C
C	[DH.3.2] COMPUTE STD DEV
	FA=SQRT(CURSYM)
C
C
C	[DH.4] NOW COMPUTE 2-SIDED STD DEV FOR OBJECT
C	COMPUTE: VARIANCE=(1/N)*SUM(IH(MEANO)-IH(I))**2 FOR I=MING:MAXG.
	IB=0
S	DCA \IB#
C
C	COMPUTE OVER 1:256 SO NO PROBLEM W/DO LOOP.
	I11=MING+1
	I12=MAXG+1
	I14=IH(MEANO+1)
	DO 410 I13=I11,I12
	IZ=I13-1
S	TAD PIH
S	CPAGE 5
S	TAD \IZ
S	DCA 7
S	TADI 7
S	CIA
S	TAD \I14 /IH(MEANO)
S	SPA /FORCE > 0
S	CIA
S	DCA \IA /CONTAINS:   ABS(IH(MEANO)-IH(I13-1)).
C
C	COMPUTE: IA=IA**2=IA*IA
S	CPAGE 12
S	TAD \IA
S	SWAB /WITH MQL
S	MUY
S	\IA /CONTAINS (IH(MEANO)-IH(I13-1))
S	DAD  /NOW COMPUTE: IB[1:2]=IB[1:2] + (IH(MEANO)-IH(I13-1))**2
S	\IB
S	DST
S	\IB
S	CLAMQ
410	CONTINUE
S	SWBA
	N=MAXG-MING+1
C
C
C	[DH.4.1] VARIANCE=IB[1:2]/N
	CALL DPCVRT(IB,CURSYM,-1)
	CURSYM=CURSYM/FLOAT(N)
C
C
C	[DH.4.2] COMPUTE STD DEV
	FB=SQRT(CURSYM)
C
C
C
C	[5] RESTORE X AND Y
500	IX=JXSAVE
	IY=JYSAVE
S	JMP RDOHIST
C	**********************************************
C	*SUBROUTINE   P R I N T H S W I T C H
C	***************************************************
S	CPAGE 3
S RPRINTHSWITCH,	JMP I PRINTHSWITCH
S PRINTHSWITCH,	0
C
C	IF /H THEN PRINT THE RESULTS
S	TAD \ISWH
S	SNA CLA
S	JMP RPRINTHSWITCH /NO, RETURN
C
	DO 1332 INDEX=1,LCNT,2
	WRITE(INDEX,332)MING,MAXG,MODEMEDIAN,MODE,FA,MEANO
     1,FB,IP,MODEDPTILE,IOPTOP,ITSLICE
332	FORMAT(' H-RNG[',I3,':',I3,'], MED=',I3,', MODE=',I3
     1,', S.D-B=',F4.1,', MEAN-O=',I3,', S.D-O=',F4.1,/,
     2,'   %TILE=',I3,', DEL%TILE=',I2,', [2.3.',I1,'], TH=',I3)
C
	WRITE(INDEX,361)JXSAVE,JYSAVE,
     1I13,I12,I11,I14,I18,I10(1),I15,I16,I17
361	FORMAT(' [',I3,',',I3,']',3(/,'     ',3I1))
C
C
C	IF /B
C		THEN PRINT BINARY VERSION OF ENTIRE WINDOW
S	TAD \ISW#
S	SNA CLA
S	JMP \1332 /NO
C
C	YES, READ IN LINE BY LINE, CVT TO 0:1 AND PRINT LINE
	DO 1333 IY1=KY1,KY2
	IY=IY1-1
	MEM=IBM1
	IBYTE=IHGH1
	CALL T3BUF(IBUF1,2)
	DO 1334 IX1=KX1,KX2
S	CLA CMA
S	TAD PBUF1
S	CPAGE 5
S	TAD \IX1
S	DCA 7
S	TADI 7
S	DCA \IZ
C	IF IZ < ITSLICE
C		THEN IZ_0 ELSE IZ_1;
S	TAD \ITSLICE
S	CIA
S	TAD \IZ
S	SMA CLA
S	IAC
S	DCA \IZ
S	CPAGE 4
S	TAD \IZ
S	DCAI 7
1334	CONTINUE
C
C	IF IY=JYSAVE
C		THEN IBUF1(JXSAVE+1)<==2+IBUF1(JXSAVE+1)
S	TAD \IY
S	CIA
S	TAD \JYSAVE
S	SZA CLA
S	JMP \1333 /NO
C
C	MARK IT
S	TAD PBUF1
S	TAD \JXSAVE
S	DCA 7
S	CPAGE 7
S	DCA \IZ /FORCE COMMON
S	TADI 7
S	TAD (2
S	DCAI 7
C
1333	WRITE(INDEX,1335)(IBUF1(IZ),IZ=KX1,KX2)
1335	FORMAT('     ',72I1)
1332	CONTINUE
C
C	RESTORE IX,IY
	IX=JXSAVE
	IY=JYSAVE
S	JMP RPRINTHSWITCH
C************** P A R A M E T E R S *************
S PI10,	\I10
S PBUF1, \IBUF1
S PIH,	\IH
	END