pr_images: color.pro

;+ PURPOSE:
;  This makes a 3-color composite image from a set of fits files.
;
; TWEAKABLE PARAMETERS:
;  the variables wr, wg, and wb adjust the color balance. They specify
;  the red, green, and blue weights.
;  rscale, gscale, and bscale determine how much the arcsinh transform
;  compresses the dynamic range. Larger values increase dynamic
;  range.
;
; DEPENDENCIES:
;  This procedure requires the IDL astronomy user's library,
;  and the David Fanning Coyote Library. 
;-
pro color

  DOSMALL = 0                   ;- use heavily binned images to get a quick preview
  MAKESMALL = 0                 ;- create the heavily binned images used by DOSMALL

  if DOSMALL then goto, small
  
;- read in images
  g = mrdfits('i.merge.smooth.fits.flat',0,h, /silent)
  r = mrdfits('r.merge.smooth.fits.flat',0,h, /silent)
  b = mrdfits('g.merge.smooth.fits.flat',0,h,/ silent)
  
;- crop to central 6K by 6K pixels
  r = float(r[1500:7499, 1500:7499])
  g = float(g[1500:7499, 1500:7499])
  b = float(b[1500:7499, 1500:7499])
  
;- make binned images, save, exit
  if MAKESMALL then begin
     r = rebin(r, 600, 600)
     g = rebin(g, 600, 600)
     b = rebin(b, 600, 600)
     save, g, r, b, file='small.sav'
     return
  endif
  
small:
  if DOSMALL then begin
     if ~file_test('small.sav') then begin
        print, 'cannot find small.sav. Run MAKESMALL=1 first.'
        return
   endif
     restore, 'small.sav'
  endif
  
;- color balance weights
  wr = .7
  wg = .7
  wb = .9

;- black threshhold
  roff = -10
  goff = -10
  boff = -10
  
;- dynamic range parameters
  rscale = 25D
  gscale = 25D
  bscale = 25D
  
;- flag bad pixels
  rbad = r eq 0
  gbad = g eq 0
  bbad = b eq 0
  
;- fill in missing regions that are present in other filters.
  lum = (r * (1 - rbad) + g * (1 - gbad) + b * (1 - bbad)) / (((1-rbad) + (1-gbad) + (1-bbad)) > 1)
  r = r * (1 - rbad) + lum * rbad
  g = g * (1 - gbad) + lum * gbad
  b = b * (1 - bbad) + lum * bbad
  
  
  bad = (r eq 0) and (g eq 0) and (b eq 0)
  bad = where(bad, ct, complement = good)
  if ct ne 0 then begin
     r[bad] = min(r[good])
     g[bad] = min(g[good])
     b[bad] = min(b[good])
  endif
  
;- arcsinh transform
  r = asinh(((r > roff) - roff) / rscale)
  g = asinh(((g > goff) - goff) / gscale)
  b = asinh(((b > boff) - boff) / bscale)
  
color:
  
  sz = size(r)
  sz = long(sz)
  
  color = fltarr(3, sz[1], sz[2])
  
;- make the color composite
  color[0,*,*] = wr * r 
  color[1,*,*] = wg * g
  color[2,*,*] = wb * b
  
;-look at the image. Requires the David Fanning COYOTE library
  tvimage, bytscl(color), /true, /keep
  
;- save images
  file='color'
  if (DOSMALL) then file+='.small'
  verbiage, 'file saving', 0, 0
  write_png, file+'.png', bytscl(color)  
end