pro ipsim10,ratio,min,max,bad_pix_map
; 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(364:*,*),2),19,'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)
p1=rotate(p1,2)
image_scale_auto,min,max,bad_pix_map,p1,pscaled,pnoisy
pnoisy=rotate(pnoisy,2)
tvfit,rebin(pscaled,2*xsize,2*ysize,/sample),19,'Before bad pixel applied'
tv,bytscl(rebin(pscaled,2*xsize,2*ysize,/sample),min=2800,max=3700,top=!d.n_colors)
tvfit,rebin(pnoisy,2*xsize,2*ysize,/sample),17,'Before bad pixel removal'
tv,bytscl(rebin(rotate(pnoisy,2),2*xsize,2*ysize,/sample),min=2800,max=3700,top=!d.n_colors)
openw,1,'landsat_with_noise' & writeu,1,pscaled & close,1
openw,1,'landsat_with_bad_pix' & writeu,1,pnoisy & close,1 & stop
bad_pix,bad_pix_map,pnoisy,p2
;p2=pnoisy
tvfit,rebin(p2,2*xsize,2*ysize,/sample),18,'After bad pixel removal'
tv,bytscl(rebin(rotate(p2,2),2*xsize,2*ysize,/sample),min=2800,max=3700,top=!d.n_colors)
p1=rotate(p2,2)
;wdelete,0
;minp1=min(p1) & maxp1=max(p1)
minp1=min(pscaled) & maxp1=max(pscaled)
;opt_stretch,p1,pout
pout=p1
tvfit,rebin(pout,1*xsize,1*ysize,/sample),1,answer+' 1',0,1100-2*270
tv,bytscl(pout,min=2800,max=3900,top=!d.n_colors)
suffix=strcompress(string(ratio)+'.'+string(min)+'-'+string(max),/remove_all)
fn='landsat.ccd.'+suffix
tiff_write,fn,rotate(rebin(pout,1*xsize,1*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)
sqrt_dmg=s(1)
p1s=byte(f(p1))
tvfit,rebin(p1s,1*xsize,1*ysize,/sample),3,'Square-rooted 3',2*xsize+10,1100-2*270
tv,bytscl(p1s,min=170,max=250,top=!d.n_colors)
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,1*xsize,1*ysize,/sample),0,'Decompressed 0',4*xsize+20,1100-2*270
tv,bytscl(pc,min=170,max=250,top=!d.n_colors)
pcb=intarr(xsize,ysize)
pcb=b(pc)
tvfit,rebin(minp1>pcb<maxp1,1*xsize,1*ysize,/sample),5,'Unsquare-rooted 5',6*xsize+30,1100-2*270
tv,bytscl(pcb,min=2800,max=3900,top=!d.n_colors)
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,1*xsize,1*ysize,/sample),6,'Dark Correction 6',20,1044-2*270-20
tv,bytscl(pcbd,min=2800,max=3900,top=!d.n_colors)
make_flat_field,bad_pix_map,pf1
;pf1=rotate(pf1,2)
tvfit,rebin(pf1,1*xsize,1*ysize,/sample),7,'Flat Field 7',2*xsize+40,1044-2*270-20
tv,bytscl(pf1,min=2500,max=3400,top=!d.n_colors)
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,1*xsize,1*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,1*xsize,1*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
pfd=pf1w
normimg=pfd(5:xsize-3,1:ysize-2)
std=stdev(normimg,flat)
pfd=float(pfd)/flat
pcbdf=pcbd
pfd=rotate(pfd,2)
pcbdf=0>fix(float(pcbd)/(pfd>0.1)+0.5)<4095
;opt_stretch,pcbdf,pout
pout=pcbdf
tvfit,rebin(pout,2*xsize,2*ysize,/sample),13,'Final Image 13',40,1044-2*270-50
tv,rebin(bytscl(pout,min=2800,max=3700,top=!d.n_colors),2*xsize,2*ysize,/sample)
fn='landsat.best.'+suffix
tiff_write,fn,rotate(rebin(pout,1*xsize,1*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
pf1dc=rotate(pf1dc,2)
p1dcf=0>fix(float(p1dc)/(pf1dc>0.1)+0.5)<4095
;opt_stretch,p1dcf,pout
;pout=p1dcf
tvfit,rebin(pout,2*xsize,2*ysize,/sample),15,'Original Image 15', $
  2*xsize+60,1044-2*270-50
tv,rebin(bytscl(p1dcf,min=2800,max=3700,top=!d.n_colors),2*xsize,2*ysize,/sample)
fn='landsat.orff.'+suffix
tiff_write,fn,rotate(rebin(pout,1*xsize,1*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)
;openw,2,'image_stats',/append
;printf,2,ratio,min,max,' ',bad_pix_map,' ',sqrt_dmg,s(1),s(2)
;close,2
bestflat=fltarr(xsize,ysize)
bestflat=0>float(pcbd)*oflat/float(p1dcf>1)<4095
tvfit,rebin(bestflat,1*xsize,1*ysize,/sample),16,'Best Flat Field 16', $
  4*xsize+80,1044-2*270-50
;kill_w
read,'Enter percentage ',per
window,1,xsize=800,ysize=700
tv,rebin(bytscl(pcbdf,min=2800,max=3700,top=!d.n_colors),2*xsize,2*ysize,/sample),410,130
tv,rebin(bytscl(p1dcf,min=2800,max=3700,top=!d.n_colors),2*xsize,2*ysize,/sample),60,130
title=string(fix(ratio),per,format="('Image Reconstruction for c = ',i1,' and ',i3,'% Bad Pixel Map')")
xyouts,30,660,title,charsize=2.5,/device
xyouts,110,100,'Uncompressed              Compressed',charsize=2.5,/device
title=string(s(2),format="('(Total Noise/Shot Noise)  = ',f4.1)")
xyouts,200,30,title,charsize=2.0,/device
xyouts,450,40,'2',/device
psave=tvrd()
fn='mtg.'+suffix+'_'+bad_pix_map
openw,1,fn
writeu,1,psave
close,1
return
end