; initialize the filenames to be strings

;r='M63r'
;g='M63v'
;b='M63b'

READ,r,PROMPT='What is the name of the red (R) image?   '
READ,g,PROMPT='What is the name of the green (V) image?   '
READ,b,PROMPT='What is the name of the blue (B) image?   '

; read them in
r=READFITS(r+'.fits')
g=READFITS(g+'.fits')
b=READFITS(b+'.fits')

; align the images
dx=0
dy=0
; correl_optimize does a really good job of auto-aligning if
; you want to do that.
;correl_optimize,r,g,dx,dy,/NUMPIX
READ,dx,dy,PROMPT='How should the green image be shifted to line up with the red image? (x,y)   '
g=SHIFT(g,dx,dy)

;correl_optimize,r,b,dx,dy,/NUMPIX
READ,dx,dy,PROMPT='How should the blue image be shifted to line up with the red image? (x,y)   '
b=SHIFT(b,dx,dy)

; increase the color saturation of red and blue relative to green
rstretch=1.
bstretch=1.

; you can have the user control color saturation with this line
;READ,rstretch,bstretch,PROMPT='Enter powers for the R and B color stretch:   '
r=r*(r/g)^rstretch
b=b*(b/g)^bstretch



; find the medians
rmed=MEDIAN(r)
gmed=MEDIAN(g)
bmed=MEDIAN(b)

; clip them so we don't have too many outliers for later statistics 
b=b>0.
b=b<4*bmed
g=g>0.
g=g<4*gmed
r=r>0.
r=r<4*rmed



; now set the lower point
b=b-MEDIAN(b)>0.
g=g-MEDIAN(g)>0.
r=r-MEDIAN(r)>0.



; calculate some statistics
x=MOMENT(b,sdev=bsig)
x=MOMENT(g,sdev=gsig)
x=MOMENT(r,sdev=rsig)

; clip them again by sigma
b=b<3*bsig
g=g<3*gsig
r=r<3*rsig


bmax=MAX(b)
gmax=MAX(g)   ; these should just be N*sigma from above but I've done 
rmax=MAX(r)   ; it this way so that number isn't always hard-coded in 
              ; too many places

; "gamma correction" --- more voodoo
gamma=1.75

; you can have the user control gamma with this line
;READ,gamma,PROMPT='Enter gamma:   '
b=(b^(1/gamma))*255/bmax^(1/gamma)
g=(g^(1/gamma))*255/gmax^(1/gamma)
r=(r^(1/gamma))*255/rmax^(1/gamma)



; finally time for some rgb stuff

a=[[[r]],[[g]],[[b]]]
outfile=''
READ,outfile,PROMPT='What would you like to name the color image?   '
write_jpeg,outfile+'.jpg',a,true=3,quality=95

END