	program halftone
	implicit NONE
c
c		June, 1985
c		R. Goldstein
c		Based on a paper by Saghri, Hou, and Tescher
c		Converts a gray-scale image to half-tone by method
c		of error propagation. Output image is shrunk or
c		expanded depending on characteristics of target device.
c		Choice of target device is given.
c
c		July, 1986
c		Fixed dimension statements in many subroutines.
c		Also added the output of the halftoned image to a file.
c
	byte inbuf(512*1024)
	integer*2 buf1(512*1024),buf2(512*1024),work(512*1024)
	integer*4 size,ri,ci	!input image parameters.nrows,ncol,size
	integer*4 rd,cd		!outpu image nrow,ncol
	real*4 l,w,aspect	!length and width of output image,aspect ratio
	integer*4 itype		!input type - already halftone=1,else0
	integer*4 rout,cout	!scaled,stretched,magnified image dimens
c
c
	external masks		!invoke block data to load mask arrays
c
c		get the imput image into memory. Can be a gray scale
c		or already halftoned image. If halftoned, skip to output
c		section
	ri = 512
	ci = 512
	size = 512		!assume the largest to start
	call getimage(size,ri,ci,inbuf,itype)
	if (itype .eq.1)THEN	!halftoned image was already input
	   rout = ri
	   cout = ci
	   call i2image(buf1,inbuf,ri,ci)
	   goto 90
	endif
c
c		convert the image to i*2
c		input to this routine is inbuf, output is in buf1
	call i2image(buf1,inbuf,ri,ci)
c
c		get target device characteristics
	call gettarget(l,w,rd,cd,aspect)
c
c		preprocess. Expand, compress, scale
c		essentially make one grey-scale pixel per output dot
c		input is in buf1, output is in buf2
	call preprocess(buf1,buf2,work,ri,ci,rd,cd,l,w,aspect,
     1			rout,cout)
c
c		Now actually create the 0-1 image. 
c		input is in buf2, output is in buf1
	call half_tone(buf2,work,rout,cout,buf1)
c
c		
90	continue
c
c		output the results. Can be to adage, file, or crt
	call out(buf1,rout,cout,inbuf)
	stop ' Normal Termination'
	end
c
	subroutine getimage(size,ri,ci,buf,itype)
	implicit NONE
	integer*4 size,ri,ci,itype
	byte buf(size*size)		!dimension here is irrelevant. 
c
c		Gets the imput image. The image can be of standard
c		type, in which case IMLOD is called, 
c		or, it can be a halftone type. Halftone images
c		are handled here.
c
	character*1 ctype
c
10	continue
	write(6,101)
101	format(' Input type of image:'/
     1		'      I	Standard image type'/
     2		'      H	Halftone image'/
     3		' Enter option (I,H) : ',$)
	read (5,102) ctype
102	format(a)
	if (ctype .eq. 'i')ctype = 'I'
	if (ctype .eq. 'h')ctype = 'H'
	if (ctype .ne.'I' .and. ctype .ne. 'H') goto 10
c
	if (ctype .eq. 'I') then
	   call imlod_2(buf,size)
	   ri = size
	   ci = size
	   itype = 0
	else if (ctype .eq.'H') then
	   call half_file_in(buf,ri,ci,size)
	   itype = 1
	else
	   write(6,104)
104	   format(' Bad input')
	   goto 10
	endif
c
	return
	end
c
	subroutine gettarget(l,w,rd,cd,aspect)
	implicit NONE
	integer*4 rd,cd
	real*4 l,w,aspect
c
c		Gives the user a choice of targets. The data values
c		are then extracted from a data statement that are
c		appropriate for the chosen target. The data items of
c		interest are:
c			l	length of target (inches or cm)
c			w	width
c			rd	number of rows
c			cd	number of columns
c			aspect	aspect ratio. Calculated
c
	common/devtype/idev,iup,jup,rowdim,coldim,mask,dthresh,
     1		onval,offval,cfact
	integer*4 idev,iup,jup,rowdim(10),coldim(10),dthresh(10),
     1		onval(10),offval(10)
	real*4 mask(4,4,10),cfact(10)
c
	real*4 dimens(4,10)	!dimensions of various devices
	integer*4 number,iopt,len,cdin
	real*4 fcd,frd
	data dimens /		!4 items/device len,wid,cols,rows
     1		9.2,6.3,80,24,	!alph crt (vt220)
     2		11.0,11.0,512,512,	!full adage
     5		14.0,11.0,132,65,	!la120, 10 pitch, wide paper
     6		8.0,5.0,799,249,	!VT240 using sixels
     7		8.5,11.0,128,87,	!la120, 8lpi, 16cpi
     8		8.0,8.0,1000,500,	!la50 - sixels
     9		16*0/
	data number/6/		!number of devices implemented
c
10	continue
	write(6,101)
101	format(' Choose a target device: '/
     1		'      1 - VT220 alpha crt'/
     2		'      2 -  adage  '/
     5		'      3 - la120,10pitch,wide paper'/
     6		'      4 - VT240 in sixels '/
     7		'      5 - la120, 8lpi, 16cpi'/
     8		'      6 - la50 - sixels '/
     9          ' Input option (1-6): ',$)
c
	read(5,102)iopt
102	format(i)
	if ( iopt .lt. 1  .or. iopt.gt.number )goto 10
c
	l = dimens(1,iopt)
	w = dimens(2,iopt)
	fcd = dimens(3,iopt)
	frd = dimens(4,iopt)
	cd = fcd
	rd = frd
c
c		override default dimensions
	if (iopt.eq.2 .or. iopt.eq.4 .or. iopt.eq.6) THEN
	   write(6,105)dimens(3,iopt)
105	   format(' Default # columns is: ',f6.0/
     1		' You can override this now. Enter new col dimension: ',$)
	   read(5,106)len,cdin
106	   format(q,i)
	   if (len .ne. 0) THEN
	      cd = cdin
	      l = l*cd/fcd
	      w = w*cd/fcd	!reduce w and rd by same percentage
	      frd = frd*cd/fcd
	      fcd = cd
	      rd = frd
	   endif
	endif
	   
	aspect = fcd*w/(l*frd)
c
	idev = iopt		!device number
	iup = rowdim(idev)	!dimension of error matrix
	jup = coldim(idev)	!ditto
c
	return
	end
c
	subroutine preprocess(buf1,buf2,work,ri,ci,rd,cd,l,w,aspect,
     1		rout,cout)
	implicit NONE
	integer*4 ri,ci,rd,cd,rout,cout
	real*4 l,w,aspect
	integer*2 buf1(ci,ri),buf2(cd*rd),work(cd*rd)
c		dimensions of buf2,work are irrelevant
c
c			Expand, compress, scale the image
c			so that one pixel corresponds to one dot
c			of output.
c			buf1 is input image, output image is returned
c			in buf2, with dimensions rout,cout
c
	integer*4 routa,couta	!intermediate dimensions
	real*4 kci,factor,fri,frd,fci,fcd
c
	if (aspect .lt. 1.) stop ' preproc Not implem. for aspect lt 1'
	write(6,101)
101	format(' Preprocessing ...')
	fri = ri
	fcd = cd
	frd = rd
	fci = ci
c	write(6,103)ri,ci,rd,cd,aspect
103	format(' preproc..ri,ci,rd,cd,aspect: ',4i5,f8.2)
c
c		Special case of halftoning to same image size
	if (aspect .eq. 1  .and. ri.eq.rd .and. ci.eq.cd) then
	    rout = ri
	    cout = ci
	    call dupit(buf1,buf2,ri,ci)
	    return
	endif
c
c		See if image has to be stretched due to aspect ratio
c		of output device not equal to 1
	if (aspect.gt.1) then
	    routa = ri
	    couta = ci*aspect
	else if (aspect .lt. 1) then
	    couta = ci
	    routa = ri/aspect
	else		!aspect = 1
	    couta = ci
	    routa = ri
	endif
	call replicate(buf1,ci,ri,work,couta,routa,aspect)
c
c		Now shrink or magnify the image 
c
	kci = aspect*ci
	if (kci.gt.cd) then
	   if (ri.gt.rd) then
	      factor = max(fri/frd,kci/fcd)
	      cout = couta/factor
	      rout = routa/factor
	      call shrink(work,couta,routa,buf2,cout,rout,factor)
	   else
	      factor = kci/fcd
	      cout = couta/factor
	      rout = routa/factor
	      call shrink(work,couta,routa,buf2,cout,rout,factor)
	   endif
	else
	   if (ri.gt.rd) then
	      factor = fri/frd
	      cout = couta/factor
	      rout = routa/factor
	      call shrink(work,couta,routa,buf2,cout,rout,factor)
	   else
	      factor = min(frd/fri,fcd/kci)
	      cout = couta*factor
	      rout = routa*factor
	      call magnify(work,couta,routa,buf2,cout,rout,factor, buf1)
	   endif
	endif
c
	return
	end
c
	subroutine half_tone(buf2,egen,rd,cd,buf1)
	implicit NONE
	integer*4 rd,cd
	integer*2 buf2(cd,rd),egen(cd,rd),buf1(cd,rd)
c
c		implements the algorithm found in Saghri, Hou, and Tescher
c		for converting a grey-scale image to half-tone image
c		Uses the generated error matrix and a mask that varies
c		with device type
c			buf2 - input image - result of preprocessing buf1
c			egen - generated error image
c			buf1 - output image. overwrites original input
c
	common/devtype/idev,iup,jup,rowdim,coldim,mask,dthresh,
     1		onval,offval,cfact
	integer*4 idev,iup,jup,rowdim(10),coldim(10),dthresh(10),
     1		onval(10),offval(10)
	real*4 mask(4,4,10),cfact(10)
c
c
	integer*2 row,col,i,j,i1,j1,value,thresh,eprop,t,maxvalue
	integer*2 one,zero,contfact
	data maxvalue/255/
c
	real*4 m11,mdiag	!for tuning experiments
	logical in_line_mult	!do do matrix multiply in line.
c
c
	thresh = dthresh(idev)
	iup = coldim(idev)
	jup = rowdim(idev)
	one = onval(idev)
	zero = offval(idev)
	m11 = mask(2,2,idev)
	mdiag = mask(1,2,idev)
c
c		this contrast factor is a kludge so that printed
c		versions will be much lighter than on the screen. The
c		reason is that a black dot when printed "shmears" over
c		more than one pixel location, giving a darker appearance to
c		the image. Essentially therefore, we make an artificially 
c		large error value in the egen matrix when a zero is placed.
c
	contfact = cfact(idev)*maxvalue
c
	do row = 1,rd
	do col = 1,cd
	   egen(col,row) = 0
	enddo
	enddo
c
	write(6,103)
103	format(' Converting to half tone ...')
c
	in_line_mult = (iup.eq.2 .and. jup.eq.2)
	do row = 1,rd
	do col = 1,cd
c
c
cc		calculate propagated error at col,row
	   eprop = 0
	   if (in_line_mult) THEN
c		hard coded matrix multipy to replace above code for 2x2 case
		if (row.ne.1 .and. col.ne.1)THEN
		   eprop = egen( col-1, row-1)*m11 +
     1		   egen( col-1, row)*mdiag +
     2		   egen( col,row-1)*mdiag
		endif
	   ELSE
c
	   	do i = 1,iup
		   do j = 1,jup
		       i1 = col - i +1
		       j1 = row - j +1
		       if (i1 .gt. 0  .and. j1.gt.0) 
     1			eprop = eprop + egen(i1,j1)*mask(i,j,idev)
		   enddo
		   enddo
	   ENDIF
c
c		set output image based on value of eprop and pixel value
c		also calculate generated error image
	   t = buf2(col,row) + eprop
	   if (t .gt. thresh) then
	      value = one
	      egen(col,row) = t - maxvalue
	   else
	      value = zero
	      egen(col,row) = t + contfact
	   endif
	   buf1(col,row) = value
c
c
c	   write(6,101)col,row,value,egen(col,row),eprop,t
101	   format(' half_tone. col,row,H,egen,eprop,t: ',6i6)
	enddo
	enddo
c
	return
	end
c
	subroutine out(buf,rout,cout,outbuf)
	implicit NONE
	integer*4 rout,cout
	integer*2 buf(cout,rout)
	byte outbuf(cout,rout)
c
c		outputs the generated half-tone image onto a device
c
	integer*4 row,col,posrow,poscol,size,kk
	character*1 odev
c
	integer*2 temp
	byte btemp(2)
	equivalence (btemp,temp)
	byte line(1024)
	integer*4 iunit
	character*80 devname
c
10	continue
	write(6,101)
101	format(' Enter the output device:'/
     1		'     F - a file '/
     2		'     A - Adage - NOT IMPLEMENTED FOR 780'/
     3		'     C - An alpha CRT'/
     4		'     4 - VT240 screen sixels'/
     5		'     L - The system lineprinter'/
     6		'     5 - An la50 connected through a VT200 terminal'/
     7		' Enter your choice: ',$)
	read(5,102)odev 
102	format(a)
	if (odev .eq. 'f')odev = 'F'
	if (odev .eq. 'a')odev = 'A'
	if (odev .eq. 'c')odev = 'C'
	if (odev .eq. '4')odev = '4'
        if (odev .eq. 'l')odev = 'L'
        if (odev .eq. '5')odev = '5'
	if (odev.ne.'F' .and. odev.ne.'A' .and. odev.ne.'C' .AND. ODEV.NE.'4'
     1		.and. odev.ne.'L' .and. odev.ne.'5') goto 10
c
c
	if (odev .eq.'F') THEN
	   call half_file_out(buf,rout,cout,outbuf)
	endif
c
	if (odev .eq.'A') THEN
c	   do row = 1,rout
c	   do col = 1,cout
c	      temp = buf(col,row)
c	      outbuf(col,row) = btemp(1)	
c	   enddo
c	   enddo
cc
c	   size = rout
c	   posrow = 1
c	   poscol = 1
c	   call putimage(outbuf,rout,cout)
	   write(6,106)
106	   format(' Adage output not implemented for 780'/
     1		' Please choose another option')
	   goto 10
	endif
c
	if (odev.eq.'C') THEN
	   call open_device(devname,iunit)
	   do row = 1,rout
	      do col = 1,cout
	         line(col)=' '
	         if (buf(col,row).eq.1)line(col)='*'
	      enddo
	      write(iunit,120)(line(kk),kk=1,cout)
120	      format(x,1024a1)
	   enddo
	   close(unit=iunit)
	endif

	IF (ODEV.EQ.'4') THEN
	        call open_device(devname,iunit)
	        call vt240_ctl(iunit,'SS')	!screen sixels on
		CALL sixel_out(iunit,BUF, COUT, ROUT)
	        call vt240_ctl(iunit,'SO')	!screen sixels off
	        close(unit=iunit)
	ENDIF
c
	if (odev.eq.'L') THEN
	   call open_device(devname,iunit)
	   do row = 1,rout
	      do col = 1,cout
	         line(col)=' '
	         if (buf(col,row).eq.1)line(col)='*'
	      enddo
	      write(iunit,120)(line(kk),kk=1,cout)
	   enddo
	   close(unit=iunit)
	endif
c
	IF (ODEV.EQ.'5') THEN
	        call open_device(devname,iunit)
	        call vt240_ctl(iunit,'PS')	!printer sixels on
		CALL sixel_out(iunit,BUF, COUT, ROUT)
	        call vt240_ctl(iunit,'PO')	!printer sixels off
	        close(unit=iunit)
	ENDIF
c
	return
	end
c
	subroutine i2image(buf1,inbuf,ri,ci)
	integer*4 ri,ci
	byte inbuf(ci,ri)
	integer*2 buf1(ci,ri)
c
	integer*2 temp		!data coercion - cvt byt to i*2 images
	byte btemp(2)
	equivalence (temp,btemp)
c
	write(6,101)
101	format(' Converting image to I*2 ...')
c
c		convert the image to i*2 - use usual trick of data coercion
	do row = 1,ri
	do col = 1,ci
	   btemp(1) = inbuf(col,row)
	   buf1(col,row) = temp
	enddo
	enddo
c
	return
	end
c
	subroutine replicate(inbuf,ci,ri,outbuf,cout,rout,aspect)
	implicit NONE
	integer*2 ci,ri,cout,rout
	integer*2 inbuf(ci,ri),outbuf(cout,rout)
	real*4 aspect
c		One dimensional stretch of the image.
c		If aspect is >1 ,then stretch in x direction by aspect
c		If aspect is <1 , then stretch in y direct by 1/aspect
c
	real*4 factor,avg,total
	integer*4 ifact,ifact1,row,col,indout,take,i
	character*1 dir
c
	if (aspect .lt. 1) stop 'replicate not implemented for k<1'
c
	if (aspect .gt. 1) then 
	   dir = 'X'
	   factor = aspect
	else if (aspect .lt. 1) then
	   dir = 'Y'
	   factor = 1./aspect
	else if (aspect.eq.1) then
	   do row = 1,rout
	   do col = 1,cout
	       outbuf(col,row) = inbuf(col,row)
	   enddo
	   enddo
	   return
	endif
	ifact = factor
	ifact1 = ifact + 1
c	write(6,101)aspect,factor,ifact,ifact1,dir
101	format(' replcat aspct,fctr,ifact,ifact1,dir: ',2f8.2,2i6,x,a1)
c
	if (dir .eq. 'X') THEN
	   do row = 1,ri
	      indout = 1
	      avg = factor
	      total = 0
	      do col = 1,ci
	         take = ifact
	         if (avg.lt.factor) take = ifact1 
	         do i = 1,take
	             outbuf(indout,row) = inbuf(col,row)
	             indout = indout+1
	         enddo
	         total = total + take
	         avg = total/col
	      enddo
c	   write(6,103)row,indout,avg,total
103	   format(' replicate..row,indout,avg,tot ',2i5,2f8.2)
	   enddo
	endif
c
	return
	end
c
	subroutine shrink(inbuf,cin,rin,outbuf,cout,rout,factor)
	implicit NONE
	integer*4 cin,rin,cout,rout
	integer*2 inbuf(cin,rin),outbuf(cout,rout)
	real*4 factor
c
c		A 2 dimensional pixel averaging scheme to shrink
c		an image from cin x rin down to cout x rout
c
	integer*4 rowind,colind,takerow,takecol,ifact,ifact1,row,col
	real*4 avgrow,avgcol,totrow,totcol
	integer*4 level,i,j
c	byte work(512,512)	!only needed if call to OUT made in last line
c
	write(6,101)
101	format(' Shrinking image ...')
c
c	write(6,102)cin,rin,cout,rout,factor
102	format(' shrink ..cin,rin,cout,rout,factor: ',4i5,f8.2)
c
	ifact = factor
	ifact1 = ifact + 1
	avgrow = factor
	totrow = 0
	rowind = 1
	do row = 1,rout
	   takerow = ifact
	   if (avgrow .lt. factor) takerow = ifact1
	   totcol = 0
	   avgcol = factor
	   colind = 1
	   do col = 1,cout
	      takecol = ifact
	      if (avgcol .lt. factor) takecol = ifact1
	      level = 0
	      do i = colind,colind+takecol-1
	      do j = rowind,rowind+takerow-1
	         level = level + inbuf(i,j)
	      enddo
	      enddo
	      level = level/(takerow*takecol)
	      outbuf(col,row) = level
	      totcol = totcol + takecol
	      avgcol = totcol/col
	      colind = colind + takecol
	   enddo
	   totrow = totrow + takerow
	   avgrow = totrow/row
	   rowind = rowind + takerow
c	write(6,104)row,avgrow,takerow
104	format(' shrink ..row,avgrow,takerow: ',i4,f8.2,i6)
	enddo
c
c	call out(outbuf,rout,cout,work)
	return
	end
c
	subroutine magnify(inbuf,cin,rin,outbuf,cout,rout,factor, buf1)
	implicit NONE
	integer*2 cin,rin,cout,rout
	integer*2 inbuf(cin,rin),outbuf(cout,rout), buf1(cout,rout)
	real*4 factor

	call expand(inbuf, cin, rin, outbuf, cout, rout, factor,buf1)
c	call interpol (outbuf, cout, rout)

	return
	end
c
	subroutine expand (inbuf,ci,ri,outbuf,cout,rout,aspect, buf1)
	implicit NONE
	integer*2 ci,ri,cout,rout
	integer*2 inbuf(ci,ri),outbuf(cout,rout), buf1(cout,rout)
	real*4 aspect
c
c		Two dimensional stretch of the image by aspect.
c		The new aditional rows and colomns are set to 0
c		to be later interpolated with the other values.
c
	real*4 factor,avg,total
	integer*4 ifact,ifact1,row,col,indout,take,i
	character*1 dir
c
	if (aspect .lt. 1) stop 'expand not implemented for k<1'
c
	if (aspect .gt. 1) then 
	   dir = 'X'
	   factor = aspect
	else if (aspect.eq.1) then
	   do row = 1,rout
	   do col = 1,cout
	       outbuf(col,row) = inbuf(col,row)
	   enddo
	   enddo
	   return
	endif
	ifact = factor
	ifact1 = ifact + 1
c	write(6,101)aspect,factor,ifact,ifact1,dir
101	format(' expand aspct,fctr,ifact,ifact1,dir: ',2f8.2,2i6,x,a1)
c
	if (dir .eq. 'X') THEN
	   do row = 1,ri
	      indout = 1
	      avg = factor
	      total = 0
	      do col = 1,ci
	         take = ifact
	         if (avg.lt.factor) take = ifact1 
	         do i = 1,take
C		     if(i.eq.1) then
	             buf1(indout,row) = inbuf(col,row)
C		     else
C		     buf1(indout,row) = 0
C		     endif
	             indout = indout+1
	         enddo
	         total = total + take
	         avg = total/col
	      enddo
c	   write(6,103)row,indout,avg,total
103	   format(' expand..row,indout,avg,tot ',2i5,2f8.2)
	   enddo
	endif
c
c	expand in the y direction now
c	Note that the outer index limit is COUT, not CI, since we have
c	already expanded in the X direction
c
	   do col = 1, cout
	      indout = 1
	      avg = factor
	      total = 0
	      do row = 1, ri
	         take = ifact
	         if (avg.lt.factor) take = ifact1 
	         do i = 1,take
C		     if(i.eq.1) then
	             outbuf(col, indout) = buf1(col, row)
C		     else
C		     outbuf(col, indout) = 0
C		     endif
	             indout = indout+1
	         enddo
	         total = total + take
	         avg = total/row
	      enddo
c	   write(6,105)row,indout,avg,total
105	   format(' expand..row,indout,avg,tot ',2i5,2f8.2)
	   enddo
c
	return
	end
c
	subroutine dupit(buf1,buf2,ri,ci)
	implicit NONE
	integer*4 ri,ci
	integer*2 buf1(ci,ri),buf2(ci,ri)
c
	integer*4 col,row
	do col = 1,ci
	do row = 1,ri
	   buf2(col,row) = buf1(col,row)
	enddo
	enddo
	return
	end
c
	subroutine putimage(outbuf,rout,cout)
	implicit NONE
	integer*4 rout,cout
	byte outbuf(cout,rout)
c
c		puts buffer to adage, starting position 1,1
c
	integer*4 image4(512),row,col
c
	do row = 1,rout
	   do col = 1,cout
	      image4(col) = outbuf(col,row)
	   enddo
c	   call ikpwr(0,1,row,cout,image4)
c	   call ikbwt
	enddo
	return
	end



	Subroutine sixel_out(iunit,BUF,COUT,ROUT)
	integer*4 iunit
C Prints bit-map images to unit 6 using Sixel format
c	This can be vt240 screen or an la50 connected throuth a vt200
C The sixel character ranges from ? (000000) to ~ (111111) for
C a total of 64 (2**6) possible 1-by-6 vertical pixel matrixs (the "Sixel").
C
C These characters have been loaded into a look up table using
C Subroutine Sixel_Lookup(Lookup_Table) 
C
C Imbin is the Image to be printed, MUST have only two grey levels
C for proper display.
C Size is the image size (eg: 128-by-128).
C Ilevel is the grey level threshold.
C Each Sixel is found from the image from 1-by-6 slice, computing the
C binary to decimal number, looking up the sixel character from
C the lookup table, printing that sixel character, continuing for
C Size times, doing a line feed (special sixel character), droping
C down 6 lines in the array, repeating until done.
C 7/85 Karl Wooledge
	Integer COUT,ROUT, Level
	INTEGER*2	BUF(COUT,ROUT)
	Byte	LevelB(4)
	Integer Sixcell			!binary to decimal variable
	Equivalence (Level,Levelb)
	Character*1  Lookup_Table(64)	!sixel charcters
	Character*3  Ans
	CHARACTER*1	LINE(1024)

	DO I = 1, 1024
		LINE(I) = '0'
	ENDDO

c	TYPE*, ROUT,COUT, BUF(1,1), BUF(20,20), BUF(50,50)
c	READ(5,310)ANS

C Load the lookup table with sixel characters (? to ~)
	Call Sixel_Lookup(Lookup_Table)

	
	ILEVEL = 0		!SET THRESHOLD


C Load Image from byte array
	Do 200 K = 1, ROUT, 6	!Sixel is a 1-by-6 vertical matrix
	   Do 150 I = 1, COUT
	      Sixcell = 0		!Reset the Sixcell
	      Ipower  = -1		!Reset the binary to decimal conversion
	      Do 140 J = K, K + 5	!Sixel is a 1-by-6 vertical matrix
	       If(J .LE. ROUT) Then	!Check for power of six
	         Ipower = Ipower + 1	!powers of two (eg: 0,1,2,3,4,5)
	         Level = BUF(I,J) 	!get image grey level
	         If(Level .NE. ILevel) Then 
	            Sixcell = Sixcell + 2**Ipower ! Range from 0 to 63
		 End if
	       End if
140	      Continue
	      Sixcell = Sixcell + 1	!Look up table range 1 to 64
	      LINE(I) = Lookup_Table(Sixcell) !Char to be printed
150	   Continue
	   WRITE(iunit,100)(LINE(II),II = 1,COUT)
100	   Format('+',1024A1,$)
	   Write (iunit,160)
160	   Format('+','-',$)	!new line
200	Continue
	

C Hold for printing
c	Read (5,310)Ans
310	Format(A)

	Return
	End


	Subroutine Sixel_Lookup(Lookup_Table)
C Loads the lookup table for Sixel graphics characters
C 7/85 Karl Wooledge

	Character*1   Lookup_Table(64)
	Character*64  Lookup_Table2

	Data Lookup_Table2(1:34) /'?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`'/
	Data Lookup_Table2(35:64)/'abcdefghijklmnopqrstuvwxyz{|}~'/

	Do 100 I = 1,64
	   Lookup_Table(I) = Lookup_Table2(I:I)
100	Continue

	Return
	End

	Block Data masks
c
c		routine to load the masks for the supported devices.
c		masks are all the same now, but the structure is here
c		to modify them individually
c
	common/devtype/idev,iup,jup,rowdim,coldim,mask,dthresh,
     1		onval,offval,cfact
	integer*4 idev,iup,jup,rowdim(10),coldim(10),dthresh(10),
     1		onval(10),offval(10)
	real*4 mask(4,4,10),cfact(10)
c
	real*4 mask1(4,4),mask2(4,4),mask3(4,4),mask4(4,4),
     1		mask5(4,4),mask6(4,4),mask7(4,4),mask8(4,4)
c
	equivalence	(mask1,mask(1,1,1)),
     2			(mask2,mask(1,1,2)),
     3			(mask3,mask(1,1,3)),
     4			(mask4,mask(1,1,4)),
     5			(mask5,mask(1,1,5)),
     6			(mask6,mask(1,1,6))
c
	data rowdim/10*2/
	data coldim/10*2/
	data dthresh/10*128/
	data cfact/10*0/
c
c	data rowdim(6)/3/   !for tuning experiments
	data rowdim(6)/2/
c
c		cfact (contrast factor) is a kludge so that printed
c		versions will be much lighter than on the screen. The
c		reason is that a black dot when printed "shmears" over
c		more than one pixel location, giving a darker appearance to
c		the image. Essentially therefore, we make an artificially 
c		large error value in the egen matrix when a zero is placed.
c

	data cfact(6)/1/
	data cfact(3)/1/
	data cfact(5)/2/
c				INSENSTITIVE TO THRESH, so leave alone
c	data dthresh(6)/240/
c
c
c		IMPORTANT - filling an array with a data statement is
c		like reading it with an unindexed read. The first index
c		gets filled the fastest. Therefore, the top row is really
c		the first column of data, etc.
c
	data mask1/ 0,   0.4,  0,  0,      !VT220 alpha CRT
     1		    0.4, 0.2,  0,  0,
     2              8*0/
	data onval(1),offval(1)/1,0/
c
	data mask2/ 0,   0.4,  0,  0,	   ! adage 
     1		    0.4, 0.2,  0,  0,
     2              8*0/
	data onval(2),offval(2)/1,0/
c
	data mask3/ 0,   0.4,  0,  0,	   !la120, 10 pitch , wide paper
     1		    0.4, 0.2,  0,  0,
     2              8*0/
	data onval(3),offval(3)/0,1/
c
	data mask4/ 0,   0.4,  0,  0,	   !vt240 sixels
     1		    0.4, 0.2,  0,  0,
     2              8*0/
	data onval(4),offval(4)/1,0/
c
	data mask5/ 0,   0.4,  0,  0,	   !la120, 8lpi, 16cpi
     1		    0.4, 0.2,  0,  0,
     2              8*0/
	data onval(5),offval(5)/0,1/
c
c	data mask6/ 0,   0.6,  0,  0,	   !la50-through vt240
c     1		    0.2, 0.1,  0,  0,
c     2	            0,   0.1,  0,  0,
c     3              4*0/
	data mask6/ 0,   0.4,  0,  0,	   !la50 tuning expt
     1		    0.4, 0.2,  0,  0,
     2              8*0/
	data onval(6),offval(6)/0,1/
c
	end
c
	subroutine vt240_ctl(iunit,operation)
	integer*4 iunit
	integer*2 operation
c
c		controls the vt240 terminal
c		R. Goldstein		August, 1985
c
c
	byte esc,dcs,st
	data esc,dcs,st / 27,144,156/
c
	if (operation.eq.'SS') THEN	!screen in sixel mode
c	   type *,'SS'
	   write (iunit,101)dcs,'1','q'
101	   format('+',10a1,$)
c
	ELSE if (operation .eq. 'SO') THEN    !screen sixels OFF
c	   type *,'SO'
	   write(iunit,101)st
c
	ELSE if (operation .eq. 'PS') THEN    !printer sixels
c	   type *,'PS'
c		turn printer controller on -- pg4-37 of vt24 ref manual
	   write(iunit,101)esc,'[','5','i'
c		printer to graphics mode pg 44 of la50 manual
	   write(iunit,101)esc,'P','q'
c
	ELSE if (operation .eq. 'PO') THEN	!prionter sixels off
c	   type *,'PO'
c			turn off printer controller mode of vt240
	   write(iunit,101)esc,'[','4','i'
c
	ELSE 
	   type *, 'Bad operation in vt240_ctl'
c
	endif
	return
	end
c
	subroutine imlod_2(buf,retsize)
	integer*4 retsize
	byte buf(1,1)
c
c		Modification History:
c			December, 1984		Use header file (RBG)
c			May, 1985  Just loads, doesn't display.
c
c		R. B. Goldstein		September, 1984
c		Loads images of different sizes from disk
c
c		Adopted from imlod by Larry Arend
c		Outline:
c			Get file name,size and desired position
c			Read the file from disk
c			Load it into the adage
c
	character filename*35
	integer*4 size		!size of image 128,256,512
	integer*4 posrow,poscol	!position to put image 11,12,13,14,21,etc.
	integer*4 nblocks	!number of blocks on disk f(size)
	integer*4 indx		!indx into linear array
	integer*4 row		!loop control variable. row number
	integer*4 x,y		!adage coordinates (pixel)
	integer*2 header(256)	!one-block header file
	logical error		!flag for no header file
	integer*4 ip		!combined row-column integer
	character headfile*35	!name of header file
	integer*4 stype		!header(4) - storage type
c
c
c
c
c		get datafile name and open the data file
3	continue
	write(6,101)
101	format(' Enter first name of image file <pix:><.dat>: ',$)
	read(5,102)lenf,filename
102	format(q,a)
	if(lenf.eq.0) then
	   type *, ' No image loaded'
	   return
	endif
	inddot = index(filename,'.')
	indcol = index(filename,':')
	if (index(filename,'.').eq.0) filename=filename(1:15)//'.dat'
	if (index(filename,':').eq.0) filename='pix:'//filename
	inddat = index(filename,'.dat')
	headfile = filename(1:inddat)//'phd'
c
	open(unit=3,name=filename,form='UNFORMATTED',
     1		access='SEQUENTIAL',status='OLD',err=900)
c
c
c		get file size and position to put picture
	call open_head(error,header,headfile)
	size = 0
	ip = 11
	if (error) then	          !get size,pos from operator
		call get_size_pos(size,posrow,poscol,ip)
	else
		size = header(1)
	endif
c
	stype = header(4)
	if( (stype.eq.0) .or. (stype.eq.2) ) then  !pixel or word format?
		call pixload_2(buf, size, posrow, poscol, stype)
	else 
		write(6,110)
110	        format(' Word images cannot be mapped..no load done')
	endif
c
c
	close(unit=3)
	retsize=size
	return
c
c
900	continue
	stop ' Error in opening image file'
	end
c
	subroutine pixload_2(image,size,posrow,poscol,stype)
c
c		reads image from disk . Does not send to adage
c		only for pixel formatted images
c
	integer*4 size,posrow,poscol,nblocks,stype
	byte image(512,size)
	integer*4 image4(512)	!a line of image converted to i*4 words
	integer*4 x,y,row,indx
c
c
	nblocks=(size*size)/512
c
c		read image from disk - one block per record
	if (stype.eq.0) then
	   do i=1,nblocks
	      read(3) (image(j,i),j=1,512)
c	      type *,' read block ',i,'  image(1,i) ',image(1,i)
	   enddo
	else		!stype must be 2 - one giant read stmt
	   read(3)image	!512x512 requires differnt read - dont know why
	endif
c
	return
	end
c
	subroutine half_file_out(buf,rout,cout,work)
	implicit NONE
	integer*4 rout,cout
	integer*2 buf(cout,rout)
	byte work(cout,rout)
c
c		outputs an already halftoned image to disk.
c
	character*35 filename
	integer*4 row,col,lenf,inddot,indcol
	integer*2 temp
	byte btemp(2)
	equivalence (btemp,temp)
c
	write(6,101)
101	format(' Enter output filename <pix:><.hft>: ',$)
	read(5,102)lenf,filename
102	format(q,a)
	if(lenf.eq.0) then
	   type *, ' No image output'
	   return
	endif
	inddot = index(filename,'.')
	indcol = index(filename,':')
	if (index(filename,'.').eq.0) filename=filename(1:15)//'.hft'
	if (index(filename,':').eq.0) filename='pix:'//filename
c
	open(unit=3,name=filename,form='UNFORMATTED',
     1		access='SEQUENTIAL',status='NEW')
c
c		convert the image to byte
	   do row = 1,rout
	   do col = 1,cout
	      temp = buf(col,row)
	      work(col,row) = btemp(1)	
	   enddo
	   enddo
c
	write(3)rout,cout
	do row = 1,rout
	   write(3)(work(col,row),col=1,cout)
	enddo
c
	close(unit=3)
	return
	end
c
	subroutine half_file_in(buf,ri,ci,size)
	implicit NONE
	integer*4 ri,ci,size
	byte buf(ci,ri)
c
c		reads an already halftoned image into memory
c
	character*35 filename
	integer*4 inddot,indcol,lenf
c
	write(6,101)
101	format(' Enter input filename <pix:><.hft>: ',$)
	read(5,102)lenf,filename
102	format(q,a)
	if(lenf.eq.0) then
	   type *, ' No image loaded'
	   return
	endif
	inddot = index(filename,'.')
	indcol = index(filename,':')
	if (index(filename,'.').eq.0) filename=filename(1:15)//'.hft'
	if (index(filename,':').eq.0) filename='pix:'//filename
c
	open(unit=3,name=filename,form='UNFORMATTED',
     1		access='SEQUENTIAL',status='OLD')
c
c
	read(3)ri,ci
	write(6,301)ri,ci
301	format('  half_file_in ri,ci ',2i5)
	call half_file_2(buf,ri,ci)	!split routine so get correct dimension
	close(unit=3)
	return
	end
c
	subroutine half_file_2(buf,ri,ci)
	implicit NONE
	integer*4 ri,ci
	byte buf(ci,ri)
c
	integer*4 row,col
c	
	do row = 1,ri
	   read(3)(buf(col,row),col=1,ci)
	enddo
c
	return
	end
c
	subroutine open_device(devname,iunit)
	implicit NONE
	character*80 devname
	integer*4 iunit
c
c		prompts for filename and returns it with a unit number
c		to main program . Opens the filename, which is assumed
c		to actually be a devicename
c
	integer*4 len
c
	iunit= 7
c
	write(6,101)
101	format(/' You can physically output to a device on another port'/
     1		' or directly to your terminal. Please enter the device'/
     2		' name (with a colon). Ie tta0:, or txc1:, or sys$print:'/
     3		' or draft_printer:, etc. For output directly to your'/
     4		' terminal, please press RETURN :  ',$)
c
	read(5,102)len,devname
102	format(q,a)
	if (len .eq. 0) devname = 'sys$output:'
c
c		open the device with large record length,so
c		output of sixels, etc can be more efficient
	open(unit=iunit,name=devname,status='UNKNOWN',recl=1200)
c
	return
	end
