pro ipsim6,ratio,min,max
; variables:
; p = the original ccd image
; p1 = the subarray of p to be processed
; xmin,ymin,xmax,ymax = the coordinates of the subarray corners
; f = the forward square root lookup table
; b = the inverse square root lookup table
; p1s = the square-rooted image
; pc = the decompressed, square-rooted image
; pcb = the unsquare-rooted, decompressed image
; dark = the dark current vector from columns 6 - 9
; pcbd = the dark-corrected, unsquare-rooted, decompressed image
; pf = the raw flat-field image pixels
; pf1 = the subarray of the flat-field to be processed
; pf1c = the decompressed flat-field image
; pf1w = the word-sized decompressed flat-field image
; darkf = the dark vector of the flat field
; pfd = the dark-corrected, decompressed flat-field image
; pcbdf = the final flat-fielded, dark-corrected, unsquare-rooted,
;         decompressed image
; p1dc = the original subarray with dark correction
; pf1dc = the original flat field with dark correction
; p1dcf = the original subarray with dark correction and flat field correction
; bestflat = pcbd/p1dcf, the ideal flat field
loadct,0
;print,'Select the image input file.'
;d_wread,h,p
d_read,'landsat',h,p
tvfit,rotate(p,2),0,'Original Frame 0',0,800
answer=""
;print,'Enter the area to be processed (hri, mri, sli).'
ymin=0
ymax=255
;read,answer
answer="hri"
case answer of
  'hri': begin
    xmin=364
    xmax=523
    end
  'mri': begin
    xmin=54
    xmax=229
    end
  'sli': begin
    xmin=234
    xmax=361
    end
else: print,'Invalid area'
endcase
xsize=xmax-xmin+1
ysize=ymax-ymin+1
p1=intarr(xsize,ysize)
p1=p(xmin:xmax,ymin:ymax)
bad_pix,'index.2500',p1,p2
p1=p2
p1=rotate(p1,2)
image_scale_auto,min,max,p1,pscaled,pnoisy
;p1=pnoisy
wdelete,0
;minp1=min(p1) & maxp1=max(p1)
minp1=min(pscaled) & maxp1=max(pscaled)
opt_stretch,p1,pout
tvfit,rebin(pout,2*xsize,2*ysize,/sample),1,answer+' 1',0,1100-2*270
suffix=strcompress(string(ratio)+'.'+string(min)+'-'+string(max),/remove_all)
fn='landsat.ccd.'+suffix
tiff_write,fn,rotate(rebin(pout,2*xsize,2*ysize,/sample),7)
p_histo,p1,x,h,2,answer+' histogram',0,70
new_sqrt,p1,f,b
; Evaluate sqrt damage
p1s=f(p1)
p1new=b(p1s)
s=snr(p1,p1new)
print,'From square rooting alone snr = ',s(0),', mean/rmse = ',s(1), $
  ', and variance = ',s(2)
p1s=byte(f(p1))
tvfit,rebin(p1s,2*xsize,2*ysize,/sample),3,'Square-rooted 3',2*xsize+10,1100-2*270
p_histo,p1s,x,h,4,'Square-rooted Image Histogram',430,70
openw,1,'compress.inp'
writeu,1,p1s
close,1
cmd='/users/ldoose/TUB/darks_test/disrcomp compress.inp compress.tmp 256 '+ $
  string(xsize)+' '+string(ratio)
spawn,cmd
cmd='/users/ldoose/TUB/software/dcs compress.tmp compress.out '+ $
  string(ysize)+string(xsize)
spawn,cmd
wdelete,0
openr,1,'compress.out'
pc=bytarr(xsize,ysize)
readu,1,pc
close,1
spawn,'rm compress.inp compress.tmp compress.out tmp.dec'
; Evaluate compression damage
s=snr(p1s,pc)
print,'From compression alone snr = ',s(0),', mean/rmse = ',s(1), $
  ,', and variance = ',s(2)
tvfit,rebin(pc,2*xsize,2*ysize,/sample),0,'Decompressed 0',4*xsize+20,1100-2*270
pcb=intarr(xsize,ysize)
pcb=b(pc)
tvfit,rebin(minp1>pcb<maxp1,2*xsize,2*ysize,/sample),5,'Unsquare-rooted 5',6*xsize+30,1100-2*270
dark=fltarr(256)
dark=p(523-291,ymin:ymax)
pcbd=pcb
for i=0,xsize-1 do pcbd(i,*)=pcb(i,*)-dark
;tvfit,rebin(pcbd,2*xsize,2*ysize,/sample),6,'Dark Correction 6',20,1044-2*270-20
print,'Select a flat-field image'
;d_wread,hf,pf
d_read,'flatff',hf,pf
pf1=intarr(xsize,ysize)
pf1=pf(xmin:xmax,ymin:ymax)
bad_pix,'index.2500',pf1,pf2
pf1=pf2
pf1=rotate(pf1,2)
tvfit,rebin(pf1,2*xsize,2*ysize,/sample),7,'Flat Field 7',2*xsize+40,1044-2*270-20
p_histo,pf1,x,h,8,'Flat Field Histogram',860,70
new_sqrt,pf1,f,b
pf1min=min(pf1)
pf1max=max(pf1)
pfs=f(pf1)
pfs=byte(pfs)
openw,1,'compress.inp'
writeu,1,pfs
close,1
cmd='/users/ldoose/TUB/darks_test/disrcomp compress.inp compress.tmp 256 '+ $
  string(xsize)+' '+string(ratio)
spawn,cmd
cmd='/users/ldoose/TUB/software/dcs compress.tmp compress.out '+ $
  string(ysize)+string(xsize)
spawn,cmd
openr,1,'compress.out'
pfsc=bytarr(xsize,ysize)
readu,1,pfsc
close,1
spawn,'rm compress.inp compress.tmp compress.out tmp.dec'
pf1w=fix(b(pfsc))
tvfit,rebin(pf1w,2*xsize,2*ysize,/sample),9,'Decompressed Flat Field 9', $
  4*xsize+60,1044-2*270-20
p_histo,fix(pf1w),x,h,10,'Decompressed Flat Field Histogram',20,0
darkf=pf(10,ymin:ymax)
pfd=pf1w
for i=0,xsize-1 do pfd(i,*)=pf1w(i,*)-darkf
tvfit,rebin(pfd,2*xsize,2*ysize,/sample),11,'Dark Corrected Flat Field 11', $
  6*xsize+80,1044-2*270-20
p_histo,fix(pfd),x,h,12,'Dark-corrected Flat Field Histogram',450,0
normimg=pfd(5:xsize-3,1:ysize-2)
std=stdev(normimg,flat)
pfd=float(pfd)/flat
pcbdf=pcbd
pcbdf=0>fix(float(pcbd)/(pfd>0.1)+0.5)<4095
opt_stretch,pcbdf,pout
tvfit,rebin(pout,2*xsize,2*ysize,/sample),13,'Final Image 13',40,1044-2*270-50
fn='landsat.best.'+suffix
tiff_write,fn,rotate(rebin(pout,2*xsize,2*ysize,/sample),7)
print,stdev(pcbdf(5:xsize-3,1:ysize-2),mean),mean
p_histo,pcbdf,x,h,14,'Final Image Histogram',880,0
p1dc=p1
for i=0,xsize-1 do p1dc(i,*)=p1(i,*)-dark
pf1dc=pf1
for i=0,xsize-1 do pf1dc(i,*)=pf1(i,*)-darkf
normimg=pf1dc(5:xsize-3,1:ysize-2)
std=stdev(normimg,oflat)
pf1dc=float(pf1dc)/oflat
p1dcf=0>fix(float(p1dc)/(pf1dc>0.1)+0.5)<4095
opt_stretch,p1dcf,pout
tvfit,rebin(pout,2*xsize,2*ysize,/sample),15,'Original Image Flat-fielded 15', $
  2*xsize+60,1044-2*270-50
fn='landsat.orff.'+suffix
tiff_write,fn,rotate(rebin(pout,2*xsize,2*ysize,/sample),7)
print,stdev(p1dcf(5:xsize-3,1:ysize-2),mean),mean
s=snr(p1dcf,pcbdf)
print,'Loss of final image compared to optimum image = ',s(0),', mean/rmse = ',$
  s(1),', and variance = ',s(2)
bestflat=fltarr(xsize,ysize)
bestflat=0>float(pcbd)*oflat/float(p1dcf>1)<4095
tvfit,rebin(bestflat,2*xsize,2*ysize,/sample),16,'Best Flat Field 16', $
  4*xsize+80,1044-2*270-50
return
end