PRO bplot_image2

;  This is the master program for plotting multiple allsky images
;  on a polar geomagnetic grid.  File image.dat contains all of
;  image attributes necessary to plot multiple images.  Procedure
;  polarplot3.pro creates an 800 by 800 byte array showing the
;  north geomagnetic polar region down to 70 degrees latitude. 
;  Procedure allsky_mag1.pro creates the necessary mapping coordinates
;  required to place each image in the 800 by 800 polar array.
;
;  Rick Doe  SRI International  8/24/94
;
;  Modified 941214 by D P Steele to require no user input.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; OPEN IMAGE DATA FILE. 
print,'OPENING IMAGE DATA FILE'
str = ''
OPENR, unit, 'images.dat',/GET_LUN  
im_number = 0
offset = 35    ;number of characters to skip in image.dat
; Read date
READF, unit, str
; Skip lines in images.dat beginning with semicolons
WHILE(STRMID(str,0,1) EQ ";") DO READF,unit, str
sdate = STRTRIM(STRMID(str,offset,STRLEN(str)-offset),2)

; Read wavelength
x = read_real(unit,offset)
wave = LONG(x)
; Read Rayleighs per data number
raydat = read_real(unit,offset)
maxi = raydat * 255.0  ; Maximum image intensity
; Read UT hours, minutes and seconds.
x = read_real(unit,offset)
hh = LONG(X)
x = read_real(unit,offset)
mm = LONG(x)
x = read_real(unit,offset)
ss = LONG(x)
; Read number of images to be plotted.
x = read_real(unit,offset)
max_images = LONG(X)

; CREATE 800 BY 800 PIXEL GEOMAGNETIC POLAR GRID
print,'CREATING POLAR GRID'
polarplot3,hh,mm,ss,grid

; Declare and initialize arrays required for image warping.
space = intarr(4)
im1 = bytarr(800,800)
im3 = bytarr(300,300)
accum  = fltarr(800,800)
count  = fltarr(800,800)
space = intarr(4)
im1(*,*) = 0
im3(*,*) = 0
accum(*,*) = 0.
count(*,*) = 0.
loadf,'mask_900.out',900,900,0,0,mask_900

; Display polar grid with geomagnetic coastlines.
window,xsi=800,ysi=800
loadlut1,'good.cm'
tvscl,grid

; START OF IMAGE UNWARP (slightly time consuming)
while (im_number lt max_images) do begin
   print,'UNWARPING IMAGE NUMBER:','  ',im_number
   grab_info,unit,filename, xdim,ydim,xc,yc,rad,maxzd,rot,$
   alt, glat, glong, hl
   print,'FILE:   ',filename
   im_number = im_number + 1
   loadb,filename,xdim,ydim,0,0,im0

;  Check for standard sized images.  If not, then
;  normalize image dimension to 300 by 300 PIXELS at
;  a radius of 125.
   if (xdim eq 300) and (ydim eq 300) and (rad eq 125) then begin
     loadb,'mask_300.img',300,300,0,0,mask_300
     mask_300 = shift(mask_300,xc-149,yc-149)
     im0 = im0 and mask_300
     find_center,im0,xc0,yc0
     im0 = shift(im0,149-xc0,149-yc0)
     xc = 149 & yc = 149
     im3 = im0
     goto, flag9
   endif   

;  Check for image data written too close to the array edge.
   space(0) = yc - rad
   if space(0) le 8 then begin
     im0 = shift(im0,0,5)
     yc = yc + 5
   endif 
   space(2) = ydim - (yc + rad)
   if space(2) le 8 then begin
     im0 = shift(im0,0,-5)
     yc = yc - 5
   endif 
   space(1) = xc - rad
   if space(1) le 8 then begin
     im0 = shift(im0,5,0)
     xc = xc + 5
   endif 
   space(3) = xdim - (xc + rad)
   if space(3) le 8 then begin
     im0 = shift(im0,-5,0)
     xc = xc - 5
   endif

;  Expand/contract image, recenter, and stuff contents
;  into a 300 by 300 array.
   im2=congrid(im0,float(xdim)*(125./float(rad)),$
   float(ydim)*(125./float(rad)))
   find_center,im2,xc,yc
   im3(24:274,24:274) = im2(xc-125:xc+125,yc-125:yc+125)

;  Reset image attributes
   xdim = 300 & ydim = 300
   xc = 149 & yc = 149
   rad = 125
   im0 = im3

; Display current allsky image and redraw 800 by 800 grid
   flag9:
   window,2,xsi=300,ysi=300,title= filename
   im3(*,149) = 255
   im3(149,*) = 255
   tvscl,im3
   im3(*,*) = 0
   window,xsi=800,ysi=800,title='PACE GEOMAGNETIC COORDINATES'
   tvscl,grid or byte(accum/(count+.0001))

; Determine polar lookup coordinates for image
   allsky_mag1,im0,xdim,ydim,xc,yc,rad,maxzd,rot,$
   alt,glat,glong,hl,hh,mm,ss,x,y

; Expand image so it fills its warped space nicely.
   print,'SMOOTHING IMAGE:','  ',im_number
   im0 = rebin(im0,900,900) * mask_900
   x = rebin(x,900,900) * mask_900 
   y = rebin(y,900,900) * mask_900
   in0 = where(im0 ne 0,count1)
   if count1 ne 0 then begin
     for i = 0L,n_elements(in0)-1L do begin
     im1(x(in0(i)),y(in0(i)))=im1(x(in0(i)),y(in0(i))) or im0(in0(i))
        if x(in0(i)) ge 1 and x(in0(i)) le 798 and y(in0(i)) ge 1 and $
        y(in0(i)) le 798 then begin
           im1(x(in0(i))-1,y(in0(i)))=im1(x(in0(i))-1,y(in0(i))) or im0(in0(i))
           im1(x(in0(i))+1,y(in0(i)))=im1(x(in0(i))+1,y(in0(i))) or im0(in0(i))
           im1(x(in0(i)),y(in0(i))-1)=im1(x(in0(i)),y(in0(i))-1) or im0(in0(i))
           im1(x(in0(i)),y(in0(i))+1)=im1(x(in0(i)),y(in0(i))+1) or im0(in0(i))
        endif
     endfor
   endif
   ind1 = where(im1 ne 0)
   for k =0L,n_elements(ind1)-1L do begin
     count(ind1(k)) = count(ind1(k)) + 1.
   endfor
   accum(*,*) = accum(*,*) + float(im1(*,*))
   im1(*,*) = 0
   tvscl,grid or byte(accum/(count+.0001))
endwhile
FREE_LUN, unit
; END OF IMAGE UNWARP 

; Display composite image 
im1 = byte(accum/(count+.0001)) or grid
tvscl,im1

; Add annotations to the composite image
boxdev, hh, mm, wave, maxi, sdate, im1, bg, c_bar
savearr,'polar.img',bg
print,'900 BY 900 BYTE ARRAY VERSION OF MERGED IMAGE'
print,'WRITTEN TO: polar.img'
; giff = query_yn("Do you want a gif image?")
giff = 'Y'
if (giff eq 'Y') then WRITE_GIF,'polar.gif',bg
print,'GIF IMAGE polar.gif WRITTEN TO DISK'

; bw = query_yn("Do you want a black and white PostScript hardcopy?")
bw = 'Y'
if (bw eq 'Y') then BEGIN
   fname = 'polar.ps'
   laser_out,hh,mm,wave,maxi,sdate,im1, c_bar, fname
   prtcom = 'lpr -l -Pps '+fname
   ;spawn,prtcom
endif

; cb = query_yn("Do you want a color PostScript file with black background?")
cb = 'Y'
if (cb eq 'Y') then $
   color_out,hh,mm,wave,maxi,sdate,im1,c_bar,'polar.cps'

end