; RGTEST.PRO
;
; Originally written to show images from Polar Camera initial field test
; at Ryning Farm Observatory
; Modified 1993/1/29 by DPS to handle images from image tree structure
; on /jasper/cnsr3_data1
; Modified 16 July 1993 from show.pro
; Modified 17 July 1993 from roofglow.pro
;
; create byte-valued input parameters
mask=BYTARR(256,256)
; load various string-valued arrays needed
s=''
!ORDER=1
; set colour table 
loadct,0

; define paths to desired images
ohrpath='/jasper/cnsr3_data1/air/c0/f0/'
ohqpath='/jasper/cnsr3_data1/air/c1/f0/'
bgdpath='/jasper/cnsr3_data1/air/c1/f1/'
; make lists of available files of desired type
ohqfiles=FINDFILE(ohqpath+'g14??060.*',COUNT=nohq)
ohrfiles=FINDFILE(ohrpath+'g14??060.*',COUNT=nohr)
bgdfiles=FINDFILE(bgdpath+'g14??060.*',COUNT=nbgd)
; keep only valid images
ohqfiles=ohqfiles(5:*) & nohq=nohq-5
ohrfiles=ohrfiles(5:*) & nohr=nohr-5
bgdfiles=bgdfiles(5:*) & nbgd=nbgd-5
; ensure file names are sorted chronologically and return UT seconds also
timesort,ohqfiles,14,7,1993,nohq,ohqJsecs
timesort,ohrfiles,14,7,1993,nohr,ohrJsecs
timesort,bgdfiles,14,7,1993,nbgd,bgdJsecs
; read approximate dark frames
rdmeandk,0,60,dk0
rdmeandk,1,60,dk1

inf=0
inl=nohq-1

; read in each image, decode the header, subtract dark frame, rebin to 
; 256x256 if necessary, scale to 8 bits and display
ip=-1
cell=16
ncell=256/cell
corr=FLTARR(ncell,ncell)
slope=corr
intcpt=corr

FOR i=inf,inl DO BEGIN
    rdkimg,ohqfiles(i),hbq,imq & hq=gethd(hbq)
    rdkimg,ohrfiles(i),hbr,imr & hr=gethd(hbr)
    rdkimg,bgdfiles(i),hbb,imb & hb=gethd(hbb)

; extract start time of image from header
    itime=BYTARR(4)
    FOR itptr=0,3 DO itime(itptr)=hq.misc.tm.(itptr)
    utsecs=TOTAL(itime*[3600,60,1,.01])

; subtract approximate dark frames
    ciq=(imq-dk0)>0
    cir=(imr-dk1)>0
    cib=(imb-dk1)>0
; approximately subtract background image
    diq=(ciq-fq*cib)>0
    dir=(cir-fr*cib)>0
; do median smoothing on difference images
    si=INTARR(2*edge,edge)
    si(0,0)=ROUND(MEDIAN(diq,3))
    si(edge,0)=ROUND(MEDIAN(dir,3))
    cib=MEDIAN(cib,3)
    IF MAX(si) LE 255 THEN si=BYTE(si)

    FOR j=0,1 DO BEGIN	; first display Q branch image, then R branch
        im=si(j*edge:j*edge+255,*)

; estimate degree of correlation of images, blockwise, to look for
; structures in OH images that don't appear in background image
        IF j EQ 0 THEN ti=ciq ELSE ti=cir
        FOR k=0,ncell*ncell-1 DO BEGIN
	    k2=k/ncell
	    k1=k-ncell*k2
	    k4=cell*k2
	    k3=cell*k1
	    subi=ti(k3:k3+cell-1,k4:k4+cell-1)
	    subb=cib(k3:k3+cell-1,k4:k4+cell-1)
	    corr(k1,k2)=CORRELATE(subi,subb)
	    lsubi=REFORM(subi,cell*cell,1)
	    lsubb=REFORM(subb,cell*cell,1)
	    res=SVDFIT(lsubb,lsubi,2)
	    intcpt(k1,k2)=res(0)
	    slope(k1,k2)=res(1)
        ENDFOR
	wid=128
	bwid=16
	WINDOW,2,XSIZE=wid+bwid,YSIZE=3*wid
	bar,bwid,wid,wid,0
	bar,bwid,wid,wid,wid
	bar,bwid,wid,wid,2*wid
	TVSCL,REBIN(BYTSCL(corr,MIN=-1,MAX=1),wid,wid,/SAMPLE),0
	XYOUTS,bwid+wid,3*wid-!D.Y_CH_SIZE,'+1',ALIGNMENT=1,/DEVICE,COLOR=0
	XYOUTS,bwid+wid,2*wid,'-1',ALIGNMENT=1,/DEVICE,COLOR=0
	XYOUTS,bwid+wid,2.5*wid,'corr',ALIGNMENT=0.5,/DEVICE,COLOR=0 $
	      ,ORIENTATION=90
	TVSCL,REBIN(BYTSCL(slope,MIN=0.0,MAX=0.1),wid,wid,/SAMPLE),1
	XYOUTS,bwid+wid,2*wid-!D.Y_CH_SIZE,'0.1',ALIGNMENT=1,/DEVICE,COLOR=0
	XYOUTS,bwid+wid,wid,'0.0',ALIGNMENT=1,/DEVICE,COLOR=0
	XYOUTS,bwid+wid,1.5*wid,'slope',ALIGNMENT=0.5,/DEVICE,COLOR=0 $
	      ,ORIENTATION=90
	TVSCL,REBIN(BYTSCL(intcpt,MIN=-100,MAX=100),wid,wid,/SAMPLE),2
	XYOUTS,bwid+wid,wid-!D.Y_CH_SIZE,'100',ALIGNMENT=1,/DEVICE,COLOR=0
	XYOUTS,bwid+wid,0,'-100',ALIGNMENT=1,/DEVICE,COLOR=0
	XYOUTS,bwid+wid,0.5*wid,'intcpt',ALIGNMENT=0.5,/DEVICE,COLOR=0 $
	      ,ORIENTATION=90
	STOP

; If autoscaling is desired, find lmaxi, otherwise use topnum
        IF topnum EQ 0 THEN BEGIN
; Identify useful min and max values for byte-scaling
            ih=HISTOGRAM(im(wprad),MIN=0,BIN=1)	; histogram of exposed image
            ihc=ih	; next line computes cumulative pixel value pdf
            FOR jj=1,N_ELEMENTS(ihc)-1 DO ihc(jj)=ihc(jj)+ihc(jj-1)
            ihc=ihc/FLOAT(MAX(ihc))		; normalize to unity
            lmaxi=MAX(WHERE(ihc LE .999))	; top 0.1% of pixels saturate
        ENDIF ELSE lmaxi=topnum
	lmaxbh=0

        ip=ip+1		; index number for screen position of image
; show image orientation, one time only
        IF ip EQ 0 THEN BEGIN
            azrad=0		; plot arrows showing the azimuth
            leftend=scl(res=120.,wwidth/2.)-twidthd/2.
            maxlen=(leftend-scl(res=120.,1.5*X_charwidth)) $
                   /scl(res=120.,(COS(azrad)+SIN(azrad))*X_charheight)
            arlen=MIN([FIX(maxlen),4])*X_charheight
            tailx=REPLICATE(scl(res=120.,0.5*X_charwidth+SIN(azrad)*arlen),2)
            taily=REPLICATE(scl(res=120.,nincol*edge),2)
            headx=[tailx(0)-scl(res=120.,arlen*SIN(azrad)), $
                   tailx(1)+scl(res=120.,arlen*COS(azrad))]
            heady=[scl(res=120.,nincol*edge+arlen*COS(azrad)), $
                   scl(res=120.,nincol*edge+arlen*SIN(azrad))]
            arrow,tailx,taily,headx,heady,/device,hsize=-0.2
            XYOUTS,headx(0)+scl(res=120.,X_charwidth) $
                  ,MIN([heady(0),scl(res=120.,wheight-X_charheight)]) $
                  ,'N',ALIGNMENT=0,/DEVICE,CHARSIZE=csz
            XYOUTS,headx(1)+scl(res=120.,X_charwidth),heady(1),'E' $
                  ,ALIGNMENT=0.5,/DEVICE,CHARSIZE=csz
        ENDIF
; display the image
        mtv,ROTATE(im*mask,5),low=lmaxbh,high=lmaxi,res=120. $
           ,(ip MOD ninrow)*edge,(nincol-1-(ip MOD ninwin)/ninrow)*edge
;   mtv,im*mask,low=lmaxbh,high=lmaxi,res=120. $
;      ,(ip MOD ninrow)*edge,(nincol-1-(ip MOD ninwin)/ninrow)*edge
; display image sequence number
        XYOUTS,scl(res=120.,(ip MOD ninrow)*edge+1) $
              ,scl(res=120.,(nincol-(ip MOD ninwin)/ninrow)*edge-1-X_charheight) $
              ,STRING(j,extension(ohqfiles(i)),FORMAT='(I1,"0:",A3)') $
	      ,/DEVICE,ALIGNMENT=0.0 $
	      ,CHARSIZE=csz,COLOR=0.8*!D.N_COLORS
; and image time
        XYOUTS,scl(res=120.,(ip MOD ninrow)*edge+1) $
              ,scl(res=120.,(nincol-1-(ip MOD ninwin)/ninrow)*edge+3+X_charheight) $
              ,/DEVICE,ALIGNMENT=0.0,CHARSIZE=csz,COLOR=0.8*!D.N_COLORS $
              ,STRING(itime(0:2),FORMAT='(I2.2,":",I2.2,":",I2.2)')
; and exposure time
        XYOUTS,scl(res=120.,(ip MOD ninrow)*edge+1) $
              ,scl(res=120.,(nincol-1-(ip MOD ninwin)/ninrow)*edge+1) $
              ,STRING(hr.misc.exp_scale*hr.exp.exposure/1000. $
              ,FORMAT='(F5.1,1X,"s")'),/DEVICE,ALIGNMENT=0.0,CHARSIZE=csz $
              ,COLOR=0.8*!D.N_COLORS
; and maximum displayed DN
        XYOUTS,scl(res=120.,((ip MOD ninrow)+1)*edge-1) $
              ,scl(res=120.,(nincol-(ip MOD ninwin)/ninrow)*edge-1-X_charheight) $
              ,STRING(FIX(lmaxi),FORMAT='("< ",I0.0," DN")') $
              ,/DEVICE,ALIGNMENT=1.0,CHARSIZE=csz,COLOR=0.8*!D.N_COLORS
; and minimum displayed DN
        XYOUTS,scl(res=120.,((ip MOD ninrow)+1)*edge-1) $
              ,scl(res=120.,(nincol-1-(ip MOD ninwin)/ninrow)*edge+1) $
              ,STRING(FIX(lmaxbh),FORMAT='("> ",I0.0," DN")'),/DEVICE $
              ,ALIGNMENT=1.0,CHARSIZE=csz,COLOR=0.8*!D.N_COLORS
        PRINT,'.',FORMAT='(A,$)'
        IF (ip+1) MOD ninwin EQ 0 THEN PRINT,'|',FORMAT='(A,$)'
        IF GET_KBRD(1) EQ 'q' THEN GOTO,DONE
    ENDFOR
ENDFOR
DONE:
PRINT,''
END