;+
; PURPOSE:
;  This procedure tries to remove additive offsets between different
;  tiles in a mosaiced GPC1 image. This seems to do a pretty good job
;  at removing such offsets in images created from the PR_DARKTEST
;  reduction class in the IPP.
;
; CATEGORY:
;  PS PR
;
; CALLING SEQUENCE:
;  tile, [file, num]
;
; INPUTS:
;  file: The filename of a mosaiced image to correct.
;  num: The number of tiles along each dimension of the mosaic. That
;       is, the mosaic is treated as a square grid of num x num
;       tiles. 
;
; OUTPUTS:
;  The procedure adds a constant to each tile in an attempt to remove
;  any offsets. It saves the result to a file named file+'.flat'.
;
; METHOD:
;  The procedure performs a least-squares minimization of the
;  discrepancy between pixel brightnesses along adjacent edges. 
;
; TWEAKABLE PARAMETERS:
;  The easiest variable to tweak is the navg variable, which changes
;  the width of the strip along each edge to use for
;  comparison. Values of 10-30 seem to work well.
;-
pro tile, file, num
  
  navg = 20                     ;-pixel width of strips to compute
  bad = 0                       ;-pixel value denoting bad or missing data

  if n_elements(file) eq 0 then file = '/media/cave/lulin/r.merge.smooth.fits'
  if n_elements(num) eq 0 then num  = 12 ;-number of tiles along each direction

  image = mrdfits(file, 0, h)
  
  sz = size(image)
  delt = sz[1]  / num           ;-pixel length of each tile
  
  wherebad = where(image eq bad)
  image[wherebad] = !values.f_nan

  tiles = fltarr(num, num)

  ;- save the edges for each tile
  edges = dblarr(4 * num^2, delt) * !values.d_nan
  for i = 0, num - 1, 1 do begin
     for j = 0, num - 1, 1 do begin
        x1 = i * delt
        x2 = i * delt + navg - 1
        x3 = (i + 1) * delt - navg
        x4 = (i + 1) * delt - 1

        y1 = j * delt
        y2 = j * delt + navg - 1
        y3 = (j + 1) * delt - navg
        y4 = (j + 1) * delt - 1

        m1 = image[x1 : x4, y1 : y2]
        m2 = image[x3 : x4, y1 : y4]
        m3 = image[x1 : x4, y3 : y4]
        m4 = image[x1 : x2, y1 : y4]
        
        ;-Option to save the entire edge profile. In practice, this 
        ;-doesn't seem to work as well as just saving the median
        ;-brightness along each edge
        ;m1 = median(m1, dimen = 2)
        ;m2 = median(m2, dimen = 1)
        ;m3 = median(m3, dimen = 2)
        ;m4 = median(m4, dimen = 1)

        m1 = median(m1)
        m2 = median(m2)
        m3 = median(m3)
        m4 = median(m4)

        edges[(i + j * num) * 4 + 0, *] = m1
        edges[(i + j * num) * 4 + 1, *] = m2
        edges[(i + j * num) * 4 + 2, *] = m3
        edges[(i + j * num) * 4 + 3, *] = m4
     endfor
  endfor
  
  a = dblarr(num^2, num^2)
  b = dblarr(num^2)
  
  ;- set up the least squares array.
  for i = 0, num - 1, 1 do begin
     for j = 0, num - 1, 1 do begin

        ij = i + j * num
        up = ij + num
        dn = ij - num
        lf = ij - 1
        ri = ij + 1
        
        ;- is there a bottom edge?
        if i ne 0 then begin
           diff = edges[lf * 4 + 1, *] - edges[ij * 4 + 3, *]
           if total(finite(diff)) ne 0 then begin
              a[ij, ij] += 1
              a[lf, ij] -= 1
              b[ij] += total(diff, /nan) / total(finite(diff))
           endif
        endif

        ;- is there a top edge?
        if i ne (num - 1) then begin
           diff = edges[ri * 4 + 3, *] - edges[ij * 4 + 1, *]
           if total(finite(diff)) ne 0 then begin
              a[ij, ij] += 1
              a[ri, ij] -= 1
              b[ij] += total(diff,/nan) / total(finite(diff))
           endif
        endif

        ;- is there a left edge?
        if j ne 0 then begin
           diff = edges[dn * 4 + 2, *] - edges[ij * 4 + 0, *]
           if total(finite(diff)) ne 0 then begin
              a[ij, ij] += 1
              a[dn, ij] -= 1
              b[ij] += total(diff, /nan) / total(finite(diff))
           endif
        endif

        ;- is there a right edge?
        if j ne (num - 1) then begin
           diff = edges[up * 4 + 0, *] - edges[ij * 4 + 2, *]
           if total(finite(diff)) ne 0 then begin
              a[ij, ij] +=1
              a[up, ij] -= 1
              b[ij] += total(diff, /nan) / total(finite(diff))
           endif
        endif

     endfor
  endfor 

  ;- solve a * x = b
  sol = la_least_squares(a, b, /double, $
                         stat = stat, $
                         residual = r)
  print, stat,r , format = '("Status: ", i0, " Residual: ", f0)'
  

  ;- see how well this flattens the 'tiles' variable
  sol = reform(sol, num, num)
  image[wherebad] = bad
  small = rebin(image, 900, 900)
  fix = rebin(sol, 900, 900, /sample)
  fix -= mean(fix)
  flattened = small + fix
  flattened[where(small eq bad)] = bad
  tvimage, bytscl(sigrange(flattened)), /keep
  
  ;- apply the correction and save
  sol -= mean(sol)
  for i = 0, num - 1, 1 do begin
     for j = 0, num - 1, 1 do begin
        x1 = i * delt
        x2 = i * delt + delt - 1
        y1 = j * delt
        y2 = j * delt + delt - 1
        image[x1 :x2, y1 : y2] += sol[i,j]
     endfor
  endfor
  ;- replace the BAD pixels
  image[wherebad] = bad
  writefits, file+'.flat', image
  
end