Pro heliumdfe
  syt0 = Systime(/seconds)

  If ( !version.os eq 'Win32' ) Then st = '\'
  If ( !version.os eq 'linux' ) Then st = '/'

  Print, 'Where is the data?'

  dir_read = DIALOG_PICKFILE( /DIRECTORY )

  ;file_mkdir,dir_read+'dfe'

  dir_write = 'J:\20120613\'
  ;dir_write = dir_read+'dfe'+st   ;First should create a directory named dfe

  ;list = FILE_SEARCH( dir_read + 'dark*.fts', COUNT = n )

  ;dark = readfits(list[0])
  ;dark = avg(dark,2)

  ;list = FILE_SEARCH( dir_read + 'flat*.fts', COUNT = n )
  ;flat = readfits(list[0])
  ;flat = avg(flat,2)
  ;flat = flat-dark

  dark = READFITS( dir_read + 'dark_0.fts' )
  flat = READFITS( dir_read + 'flat_0.fts' )

  sm = 5.
  ns = 3.

  dirt = flat - SMOOTH( flat, sm, /EDGE )
  sigma = STDEV(dirt, mean)
  dirt1 = (dirt GE mean + ns * sigma) or (dirt LE mean - ns * sigma)
  dirt = DILATE(dirt1,REPLICATE(1,2,2))
  dirt_ind = WHERE( dirt EQ 1 )



  index1 = WHERE( flat LE 0.0, count )
  index2 = WHERE( flat GT 0.0 )
  IF ( count GT 0 ) THEN flat( index1 ) = avg( flat( index2 ) )

  ;flat = flat/avg(flat)

  list = FILE_SEARCH( dir_read + 'i*.*', COUNT = nn )

  ;IF NOT KEYWOD_SET(sm) THEN sm = 10
  ;IF NOT KEYWOD_SET(ns) THEN ns = 5


  tmp_img = READFITS( list[0] )
  s = Size(tmp_img,/Dimensions)

  xx = s[0]
  yy = s[1]
  zz = s[2]

;Following create the circle index of the bad pixels
  xt=transpose(replicate(1,yy)) ## indgen(xx)
  yt=indgen(yy) ## replicate(1,xx)
  xcen=801
  ycen=263
  rad=12
  indc=where( sqrt((xt-xcen)^2+(yt-ycen)^2) lt rad)

  sample_img = FLOAT( tmp_img[ *, *, 0 ] - dark ) / FLOAT( flat )

  Xtvscl, sample_img >0<6000

  BOX_CURSOR, xo, yo, deltx, delty, /MESSAGE

  con = FLTARR(zz)

  spos = Strpos( list[0], st, /reverse_search )
  sspos = Strpos( list[0],'_', /reverse_search )
  dpos = Strpos( list[0], '.' )
  sn = dpos - spos
  snn = dpos - sspos
  img_cube = FLTARR(xx, yy, zz)

FOR i = 0, nn-1 DO BEGIN
  ;FOR i = 0, nn - 1 DO BEGIN

    stime0 = Systime(/seconds)

    img = READFITS( list[i], hdr )

    Print,i

    FOR j = 0, zz - 1 DO BEGIN

      imgg = FLOAT( img[ *, *, j ] - dark ) / FLOAT( flat )
      area_s = imgg[ xo : xo + deltx, yo:yo + delty]

      con(j)=stdev( area_s )/Avg( area_s )

      ;index1 = WHERE( imgg LE 0.0 or imgg gt 8000., count )
      ;index2 = WHERE( imgg le 8000. and imgg gt 0.0)
      ;IF ( count GT 0 ) THEN imgg( index1 ) = avg( imgg( index2 ) )

     med = smooth(imgg,40,/edge)
      imgg[indc] = med[indc]

 IF ( N_ELEMENTS( dirt_ind )  NE 1 ) THEN BEGIN

         med = MEDIAN(imgg, 5)
        imgg(dirt_ind) = med(dirt_ind)

 ENDIF


      imggg = cao_erasedirt(imgg, 5, 1., 3.)

 IF ( N_ELEMENTS( dirt_ind )  NE 1 ) THEN BEGIN

         med = MEDIAN(imggg, 5)
        imggg(dirt_ind) = med(dirt_ind)

 ENDIF
;Correct the bad column
      smimg = smooth(imggg[253:264,*],4,/edge)
      imggg[254:261,*]= smimg[1:8,*]
      smimg = smooth(imggg[509:518,*],4,/edge)
      imggg[510:515,*]= smimg[1:6,*]
      smimg = smooth(imggg[765:774,*],4,/edge)
      imggg[766:771,*]= smimg[1:6,*]
      imggg[*,708] = (imggg[*,706]+imggg[*,707]+imggg[*,709]+imggg[*,710])/4.


      img_cube(*, *, j) = imggg

      TVSCL, imggg >2000<6000



    ENDFOR

    s = REVERSE( SORT(con) )

    sel_img = img_cube(*,*,s[0])

    print,s[0]
    tvscl,img_cube(*,*,s[0]) >2000<6000


    fname = Strmid( list[i], spos + 1, sn - 1 )

    WRITEFITS, dir_write + fname + '_sel.fts', sel_img

    OPENW, 1, dir_write + 'sel.' + num2str(trim(i,1), length=3, padchar='0')
    WRITEU, 1, img_cube
    CLOSE, 1

    OPENW, 1, dir_write + 'sel.' + num2str(trim(i,1), length=3, padchar='0') +'.txt'
    PRINTF, 1, hdr
    CLOSE, 1

    stime1 = Systime(/seconds)
    dt = stime1 - stime0
    adt = Fix(dt * Float( nn - 1 - i ))

    Print, 'Image #' + Strtrim(String(i+1),2) + ' / ' + Strtrim(String(nn),2)
    Print, Strtrim( String(adt), 2 ) + ' Seconds left'


  ENDFOR

syt1 = Systime(/seconds)
dsyt = syt1 - syt0

PRINT,'Work has been done! (' + Strtrim( String( dsyt ), 2 ) + ' seconds )'



END