C	PROGRAM SHADE.FT
C	----------------
C
C	PETER LEMKIN AND MORTON SCHULTZ
C	IMAGE PROCESSING UNIT, DCBD
C	NATIONAL CANCER INSTITUTE
C	NATIONAL INSTITUTES OF HEALTH
C	9000 ROCKVILLE PIKE
C	BETHESDA, MD. 20014
C
C	APRIL 8, 1978
C	APRIL 7, 1978
C	APRIL 6, 1978
C	APRIL 4, 1978
C	FEB 7, 1978
C	JAN 27, 1978
C
C	INTRODUCTION
C	------------
C	SHADE.FT IS A CHAINED PROGRAM USED WITH BMON2
C	IT RECEIVES ITS ARGUMENTS FROM THE CD AREA AND THE IBM1,IHGH1
C	VARIABLES IN COMMON.
C
C	<BMj>_SHADE,<BMi1>,<OPT. BMi2> - shade correct BMi into BMj given
C	a set of Neutral density images in BMi2 (or BM3 if BM2
C	is not specified). The algorithm was specified by Mort Schultz
C	as requiring the following setup procedure.
C
C	[1] All operations must be done w/QMT "light sens" switch in
C	the manual position. Adjust the image on the TV system using
C	the QMT "sensitivity" control to get an image
C	with a good dynamic range. The meter should read from 0.3
C	to 0.4.
C
C	[2] Position all BMs in the same F&S position to acquire data
C	using POSFS or ALL384.
C
C	[3] Acquire the two shading ND filter imges with no
C	specimen in the system. (Clear field)
C
C		[3.1] Insert ND filter #1 (dark filter) and get
C		data into the low BM as:   GET,BM3.
C
C		[3.2] Insert ND filter #2 (light filter) and get
C		data into the high BM as:    GET,BM3H.
C
C	[4] Without moving either the QMT sensitivity control the
C	BM positions (relative to the LCS camera position) acquire one or
C	more data images at the same position on the camera. This implies
C	that either the stage or film be moved for the microscope
C	or film reader respectively.
C
C	[5] Run the SHADE operator to shade correct the data image.
C
C	It first computes Iavg1 and Iavg2 from the Logs of two ND images 
C	where image 2 is lighter than image 1.
C
C	Let FC=255.0/Log(256).
C
C	Then,
C	Iavgk =   1         * SUM (A*Log(Gndk(x,y)+1)
C	  	-----------
C		A*(TOTN)  x,y in [0:255]
C
C	Then
C
C	K=Iavg1-Iavg2.
C
C	Finally,
C
C		         -1
C	Ishade(x,y) = Log  (K * (G(x,y)-Gnd1(x,y))  + Iavg1) - 1;
C			        ------------------
C			           (Gnd2-Gnd1)
C	OPDEFS
C	------
C
S	OPDEF	LDXP	6443
S	OPDEF	LDYP	6444
C
S	OPDEF	TADI	1400
S	OPDEF	DCAI	3400
C
S	OPDEF	DISP1	6435
S	OPDEF	DISP2	6436
C
C
S	OPDEF	SWBA	7447
C
S	OPDEF	MUY	7405
S	OPDEF	DVI	7407
C
S	OPDEF	MQA	7501
S	OPDEF	MQL	7421
S	OPDEF   KRS	6034
S	OPDEF	BSW	7002
C
C
C
	DIMENSION ILOG(256),IALOG(256)
C	[1] INITIALIZATION
	WRITE(1,995)
	WRITE(3,995)
995	FORMAT('0 SHADE 4/8/78 - 2:37PM')
	FC=TIMER(0)
C
C
C	[1.1] GET THE BMi BM NAME
S	TAD (4040 /"  "
S	DCA \I11
S	TAD \IHGH1
S	RAR /BIT 11==>LINK
S	CLA
S	TAD (1000 /"H@"
S	SZL /ADD "H" IF HIGH PART
S	DCA \I11
S	CLA
C
C
C	[1.2]GET THE BMj BM NAME
S	TAD (4040 /"  "
S	DCA \I12
S	TAD \JHGH
S	RAR /BIT 11==>LINK
S	CLA
S	TAD (1000 /"H@"
S	SZL /ADD "H" IF HIGH PART
S	DCA \I12
S	CLA
C
C
C	[1.3] PRINT BM NAME INFORMATION
	DO 250 INDEX=1,3,2
250	WRITE(INDEX,251)JBM,I12,CURSYM,IBM1,I11,KX1-1,KX2-1,KY1-1
     1,KY2-1
251	FORMAT(' BM',I1,A1,'<==' BM',I1,A1,' [',I3,':'
     1,I3,' , ',I3,':',I3,']')
C
C
C	[1.4] GET THE BMJ POSITION
C	FROM COMMON - SET BY BMON2.SV 
	IXPOSITION=LSAVE(13,JBM+1)
	IYPOSITION=LSAVE(14,JBM+1)
C	[2.0] VERIFY BM SPECS
S	JMS CKIN
S	JMS CKOUT
C
C
C	[2.1] SETUP THE SHADING MEMORY
C	DEFAULT TO BM3 IF NOT SPECIFIED
S	TAD \SEXT
S	CIA
S	TAD BMTEXT
S	MQA
S	SNA CLA
S	JMP \210
	IBM2=3
	GO TO 211
S \210,	JMS CKIN2
C
211	ISBM1=IBM2
	ISBY1=1
	ISBM2=ISBM1
	ISBY2=0
C
C
C	[2.2] COMPUTE IAVG1 AND IAVG2
C	ZERO DOUBLE PRECISION COUNTERS
	IA=0
	IB=0
	IC=0
S	DCA \IA#
S	DCA \IB#
S	DCA \IC#
C
C
C	PROCESS THE IMAGE
	DO 200 IY1=KY1,KY2
C
S	CPAGE 3
S	JMS TTYCTL	/TEST FOR TTY INPUTS
S	JMP \998	/ERROR RETURN: ABORT
C
	IY=IY1-1
S	TAD \IY
S	DISP2		/PUT BMY ADDRESS IN RIGHT LED'S
C
C	LOAD THE Y CURSOR
S	TAD \IY
S	TAD \IYPOSITION
S	LDYP
C
C	READ THE ND IY LINES
C	ND1 IN IBUF1
	MEM=ISBM1
	IBYTE=ISBY1
	CALL T3BUF(IBUF1,2)
C
C	ND2 IN IBUF2
	MEM=ISBM2
	IBYTE=ISBY2
	CALL T3BUF(IBUF2,2)
C
C
C	[2.2.1] PROCESS A LINE
	DO 200 IX1=KX1,KX2
S	CLA CMA
S	TAD \IX1
S	DCA \IX
C
S	TAD \IX
S	DISP1
C	LOAD THE X CURSOR
S	TAD \IX
S	TAD \IXPOSITION
S	LDXP
C
C	INCREMENT TOTAL NUMBER OF PIXELS
S	ISZ \IC
S	SKP
S	ISZ \IC#
S	NOP
C
C
C	[2.2.1.1] COMPUTE THE RUNNING SUMS OF ND1, ND2
S	TAD PBUF1
S	CPAGE 5
S	TAD \IX
S	DCA 7
S	TADI 7
S	DCA \IVAL
S	JMS LOKLOG /IVAL<==IFIX(FB*LOG(IFLOAT(IVAL+1))
S	TAD \IVAL
S	CLL
S	TAD \IA
S	DCA \IA
S	RAL
S	TAD \IA#
S	DCA \IA#
C
S	TAD PBUF2
S	CPAGE 5
S	TAD \IX
S	DCA 7
S	TADI 7
S	DCA \IVAL
S	JMS LOKLOG /IVAL<==IFIX(FB*LOG(IFLOAT(IVAL+1))
S	TAD \IVAL
S	CLL
S	TAD \IB
S	DCA \IB
S	RAL
S	TAD \IB#
S	DCA \IB#
C
C
200	CONTINUE
C
C
C	[2.3] CVT D.P. SUMS TO FLOATING POINT
	CALL DPCVRT(IA,FA,-1)
	CALL DPCVRT(IB,FB,-1)
	CALL DPCVRT(IC,FC,-1)
C
C
C	[2.3.1] COMPUTE THE AVERAGES
	IAVG1=(1.0/(FC*FC))*FA
	IAVG2=(1.0/(FC*FC))*FB
C
C
C	[2.3.2] COMPUTE THE AVG DIFFERENCE
	K=IAVG1-IAVG2
C	[3] SHADE CORRECT THE IMAGE
	DO 300 IY1=KY1,KY2
C
S	CPAGE 3
S	JMS TTYCTL	/TEST FOR TTY INPUTS
S	JMP \998	/ERROR RETURN: ABORT
C
	IY=IY1-1
S	TAD \IY
S	DISP2		/PUT BMY ADDRESS IN RIGHT LED'S
C	LOAD THE Y CURSOR
S	TAD \IY
S	TAD \IYPOSITION
S	LDYP
C
C
C	[3.1] GET THE THREE INPUT LINES
C	ALSO GET AN OUTPUT LINE IN CASE SHADE CORRECTION IS FOR LESS THAN
C	A FULL 256X256 IMAGE.
C
C	ND1 IN IBUF1
	MEM=ISBM1
	IBYTE=ISBY1
	CALL T3BUF(IBUF1,2)
C
C	ND2 IN IBUF2
	MEM=ISBM2
	IBYTE=ISBY2
	CALL T3BUF(IBUF2,2)
C
C	BMi IN IBUF3
	MEM=IBM1
	IBYTE=IHGH1
	CALL T3BUF(IBUF3,2)
C
C	BMj IN IBUF4
	MEM=JBM
	IBYTE=JHGH
	CALL T3BUF(IBUF4,2)
C
C
C	[3.1.1] PROCESS A LINE
	DO 310 IX1=KX1,KX2
S	CLA CMA
S	TAD \IX1
S	DCA \IX
C
S	TAD \IX
S	DISP1
C	LOAD THE X CURSOR
S	TAD \IX
S	TAD \IXPOSITION
S	LDXP
C
C
C	[3.1.1.1] LOOKUP (ND1,ND2,BMi) AT [X,Y].
C	COMPUTE: I11=IBUF1(IX1)
S	TAD PBUF1
S	CPAGE 5
S	TAD \IX
S	DCA 7
S	TADI 7
S	DCA \I11
C
C	COMPUTE: I12=IBUF2(IX1)
S	TAD PBUF2
S	CPAGE 5
S	TAD \IX
S	DCA 7
S	TADI 7
S	DCA \I12
C
C	COMPUTE: I13=IBUF3(IX1)
S	TAD PBUF3
S	CPAGE 5
S	TAD \IX
S	DCA 7
S	TADI 7
S	DCA \I13
C
C
C	[3.1.1.2] COMPUTE SHADE CORRECTED VALUE.
C	COMPUTE IZ=ANTILOG(K*(I13-I11)/(I11-I12)+IAVG1)-1)
C	COMPUTE: IVAL=K*(I13-I11)/(I11-I12)+IAVG1
S	TAD \K
S	DCA M1
S	TAD \I11
S	CIA
S	TAD \I13
S	MQL
S	TAD \I12
S	CIA
S	TAD \I11
S	DCA D1
S	CPAGE 4
S	MUY
S M1,	0	/(I13-I11)
S	DVI
S D1,	0	/(I11-I12)
S	CLA
S	MQA
S	TAD \IAVG1
S	DCA \IVAL
S	TAD \IVAL
S	JMS LOKALOG
	IZ=IVAL-1
C
C	COMPUTE: 	IBUF4(IX1)=IZ
S	TAD PBUF4
S	TAD \IX
S	DCA 7
S	CPAGE 4
S	TAD \IZ
S	DCAI 7
C
310	CONTINUE
C
C
C	[3.1.2] WRITE OUT LINE
	CALL T3BUF(IBUF4,3)
C
300	CONTINUE
C
C
C	[999] RETURN
	GOTO 998
999	WRITE(1,996)
996	FORMAT('BAD BM SPECIFICATION!')
998	FC=TIMER(2)
C
C
	CALL CHAIN('BMON2')
C***************************************************
C	*PROCEDURE      C K O U T
C******************************************************
C
C
C	CHECK WHETHER THE OUTPUT BM SPEC IS LEGAL ELSE GOTO 999.
S	CPAGE 3
SRCKOUT,	JMP I CKOUT
S CKOUT,	0 /ENTRY
C
C	[1] TEST IF KOUTFILE="BM"
S	TAD \KOUTFILE
S	CIA
S	TAD BMTEXT
S	SZA CLA
S	JMP \999 /FAILED
C
C	[2] TEST IF (KOUTFILE(2) LAND '7700)=DIGIT
S	TAD \KOUTFILE#
S	BSW
S	AND (0077
S	TAD (-60 /"0"
S	SPA
S	JMP \999 /NO, FAILED
S	TAD (-D8 /"7" TEST
S	SMA CLA
S	JMP \999 /FAILED
S	JMP RCKOUT /OK.
C
C
C***************************************************
C	*PROCEDURE      C K I N
C******************************************************
C
C
C	CHECK WHETHER THE INPUT BM SPEC IS LEGAL ELSE GOTO 999.
S	CPAGE 3
SRCKIN,	JMP I CKIN
S CKIN,	0 /ENTRY
C
C	[1] TEST IF BMI1="BM"
S	TAD \SFILE
S	CIA
S	TAD BMTEXT
S	SZA CLA
S	JMP \999 /FAILED
C
C	[2] TEST IF (BMI1(2) LAND '7700)=DIGIT
S	TAD \SFILE#
S	BSW
S	AND (0077
S	TAD (-60 /"0"
S	SPA
S	JMP \999 /NO, FAILED
S	TAD (-D8 /"7" TEST
S	SMA CLA
S	JMP \999 /FAILED
S	JMP RCKIN /OK.
C
C
C***************************************************
C	*PROCEDURE      C K I N 2
C******************************************************
C
C
C	CHECK WHETHER THE INPUT BM SPEC IS LEGAL ELSE GOTO 999.
S	CPAGE 3
SRCKIN2,	JMP I CKIN2
S CKIN2,	0 /ENTRY
C
C	[1] CHECK FIRST IBM1 SPEC
S	JMS CKIN
C
C
C	[2] TEST IF BMI2="BM"
S	TAD \SEXT
S	CIA
S	TAD BMTEXT
S	SZA CLA
S	JMP \999 /FAILED
C
C	[3] TEST IF (BMI2((6) LAND '7700)=DIGIT
S	TAD \SEXT#
S	BSW
S	AND (0077
S	TAD (-60 /"0"
S	SPA
S	JMP \999 /NO, FAILED
S	TAD (-D8 /"7" TEST
S	SMA CLA
S	JMP \999 /FAILED
S	JMP RCKIN2 /OK.
C
C
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   L O K L O G
C	********************************************************
C
C	LOOKUP FC*LOG(IVAL+1)==>IVAL
C
C	NOTE: THE RANGE OF THE INPUT IS 0:255, THE RANGE OF THE
C	OUTPUT IS 0:255.
C
S	CPAGE 3
S RLOKLOG,	JMP I LOKLOG
S LOKLOG,	0 /ENTRY
C
C	[1] IF TABLE NOT HERE THEN GENERATE IT FIRST TIME
S	TAD \NLSW
S	SZA CLA
S	JMP \1501	/NOT THE FIRST TIME
	FC=255.0/ALOG(256.0)
	DO 1500 L=1,256
	LLL=FC*ALOG(FLOAT(L))
C	MAX OF 255
S	TAD \LLL
S	TAD (-D256
S	SMA 
S	CLA CMA /-1
S	TAD (D256 /CLIP AT 255
S	DCA \LLL
1500	ILOG(L)=LLL
	NLSW=1
C
C
C	[2] DO LOOKUP
C	NOTE: VALUES OF L IN THE RANGE 0:255 COVER SUBSCRIPTS 1:256
1501	L=IVAL
C	FORCE CURRENT FIELD USING L
S	TAD \L
S	TAD PILOG
S	DCA 7
S	TADI 7
S	DCA \IVAL
C
C
S	JMP RLOKLOG /RETURN
C	********************************************************
C	*SUBROUTINE   L O K A L O G
C	********************************************************
C
C	LOOKUP IVAL<==ALOG(IVAL)
C	NOTE: THE RANGE OF THE INPUT IS 1:256, THE RANGE OF THE OUTPUT 
C	IS 0:255.
C
S	CPAGE 3
S RLOKALOG,	JMP I LOKALOG
S LOKALOG,	0 /ENTRY
C
C	[1] IF TABLE NOT HERE THEN GENERATE IT FIRST TIME
S	TAD \NALSW
S	SZA CLA
S	JMP \1601	/NOT THE FIRST TIME
	IVSV=IVAL
C
C	INIT THE TABLE
C
C	ZERO THE TABLE
	DO 1605 I13=1,256
1605	IALOG(I13)=0
C
C	COMPUTE ENTRIES WITH INTERPOLATION BETWEEN POINTS.
	DO 1602 I13=1,256
	IVAL=I13-1
S	JMS LOKLOG
S	TAD \IVAL
S	SPA /TEST IF < 0
S	JMP \1609 /YES, < 0, ERROR!
C	TEST IF > 255
S	TAD (-D256
S	SMA CLA
S	JMP \1609 /ERROR!
C
C	COVER 1:256 OF THE IALOG ARRAY  RATHER THAN 0:255
	IALOG(IVAL+1)=I13-1
C
C	SEE IF INTERPOLATE FROM LAST IVAL+1 TO CURRENT IVAL-1
	IF(I13-1)1602,1602,1603
1603	I12=IVAL-1
	DO 1604 I13=I11,I12
1604	IALOG(I13+1)=I13-1
1602	I11=IVAL+1
C
	NALSW=1
	IVAL=IVSV
C
C
C	[2] DO LOOKUP
1601	L=IVAL
C	FORCE CURRENT FIELD USING L
S	CLA CMA
S	TAD \L
S	TAD PIALOG
S	DCA 7
S	TADI 7
S	DCA \IVAL
C
C
S	JMP RLOKALOG /RETURN
C
C	ERROR! ON INITIALIZING TABLE
1609	WRITE(3,1619)L,IVAL(1)
1619	FORMAT(' ERROR OF LOKLOG(',I5,')=',I5)
	CALL CHAIN('BMON2')
C************** P A R A M E T E R S *************
S PBUF1,	\IBUF1
S PBUF2,	\IBUF2
S PBUF3,	\IBUF3
S PBUF4,	\IBUF4
S PI10,	\I10
S BMTEXT,	TEXT /BM/
S PILOG,	\ILOG
S PIALOG,	\IALOG
	END