function iris_sji_mpsf_rl1frame,index,data,psffile=fpsf,niter=niter

;+
; Name: iris_sji_mpsf_rl1frame
;
; Purpose: Deconvolve one SJI NUV frame with a spatially dependent PSF.
;          Requires a grid of PSFs stored in a named .geny file.  Each portion
;          of the image is deconvolved with the closest PSF in the grid.
;          The input image can be Level 1+ (dark subtracted, flat fielded,
;          possibly despiked and/or dustbusted) or Level 2.
;          The deconvolution uses Richardson Lucy.
;
; Ccalling sequence:
;       image_out = iris_sji_mpsf_rl1frame(index,image_in)
;
; Input:
;       index    = index of input image
;       image_in = input image
;       
; Output
;       image_out = Multi-PSF, Richardson-Lucy deconvolved image
;
; Keywords input:
;       psffile = name of .geny file with PSFs.
;                 Default: $IRIS_RESPONSE/iris_sji_nuvpsf_2016grid_l2.geny
;       niter   = number of R-L iterations.
;                 Default = ceil(2.5 * mean(data)^0.4)
;                 Note: if niter eq 0 then default is calculated and
;                       returned in niter.  Careful using niter=0 in a loop!
;                 
; Keywords output:
;
; Common blocks:
;       compsf  holds s_psf (tags xc,yc,sz,psf,sltxy),cfpsf across calls
;
; Restrictions:
;       Assumes grid of PSFs to be regular
;
; Modification history:
;       Written 2020-11-21 JPW iris_multi_psf_rl.pro   for L1 data
;               2021-02-06 JPW iris_sji_mpsf_rl1frame  for L2 and L1 data
;               2021-03-18 JPW changed restore to restgenx and
;                              added signal dependent default for Niter
;               2022-12-09 JPW Changed algorithm to find area to deconvolve,
;                              added forgiving treatment of missing pixels
;                              and not deconvolving narrow strips of image.
;                              Added lower limit of iterations in auto mode.
;               2023-01-27 JPW fixed issue with min(data) lt -200,
;                              added range check for dx1,dx2,dy1,dy2
;                              added adjustment for width lt 128 and check.
;-

common compsf,s_psf,cfpsf

maxmis = 0.5   ; maximum fraction of missing values allowed in tiles

mind = -200    ; value assigned to missing data. Needs to match value in
               ; irisl12_savesjidata.pro! Not given in header, it seems.

if n_elements(fpsf) ne 1 then begin
  respdir = get_logenv('$IRIS_RESPONSE') ; default is in $SSW_IRIS distribution
  if respdir eq '' then respdir=concat_dir('$SSW_IRIS','response')
  fpsf = concat_dir(respdir,'iris_sji_nuvpsf_2016grid_l2.geny')      ; PSF file
endif

; check if SJI NUV data
sjinuv = 0
if index.lvl_num lt 2 then begin
   if strmid(index.img_path,0,5) eq 'SJI_2' then sjinuv=1
endif else begin
   if strmid(index.tdesc1,0,5) eq 'SJI_2' then sjinuv=1
endelse
if sjinuv ne 1 then begin
   print,'Error, data is not SJI NUV image'
   return,0
endif
if index.sumspat ne 1 or index.sumsptrl ne 1 then begin
   print,'Error, not full resolution data'
   return,0
endif

; preload/restore all the PSFs from single save file
if n_elements(cfpsf) ne 1 then cfpsf=''
if cfpsf ne fpsf then begin
   restgenx,s_psf,file=fpsf
   cfpsf = fpsf
endif
sz = s_psf.sz   ; for easier reference

; replace NaNs with -200 (Note: missing L2 data usually is -200 already)
wnan = where(data ne data,nwnan)
if nwnan gt 0 then data[wnan] = mind

; make sure -200 is the lowest data value
data >= mind

; find image area to process
; first find the smallest rectangular area that includes all valid data
xok = where(max(data,dim=2) gt mind)
yok = where(max(data,dim=1) gt mind)
dx1 = min(xok)
dx2 = max(xok)
dy1 = min(yok)
dy2 = max(yok)
; second clip borders:
; Note: may no longer be needed due to tile-wise treatment of missing pixels.
;       However, removal would need to be rather thoroughly tested,
;       indluding side effects of deconvolving data with missing pixels.
; outermost 4-6 pixels may be partly invalid due to image rotation in dewarp
mar = 4+fix(((dx2-dx1)>(dy2-dy1))/280)   ; adjust margin for FOV
dx1 = dx1+mar
dx2 = dx2-mar
dy1 = dy1+mar
dy2 = dy2-mar

; If image area narrower than 128, increase if possible.
; May happen with missing data.
if dx2-dx1 lt sz-1 then begin
   ddx = sz-1-(dx2-dx1)         ; missing width
   siz_1 = (size(data))[1]-1    ; max dx2
   mrl = dx1                    ; available space on left
   mrr = siz_1-dx2              ; available space on right
   if mrr gt mrl then begin     ; more space on right
      dx1 = 0 > (dx1-ddx/2)     ; add up to ddx/2 on left
      dx2 = siz_1 < (dx1+sz-1)  ; add rest on right
   endif else begin
      dx2 = siz_1 < (dx2+ddx/2) ; add up to ddx/2 on right
      dx1 = 0 > (dx2-sz+1)      ; add rest on left
   endelse
endif

; adjust xc yc for psf slit center and data slit center
if index.lvl_num lt 2 then sltxy = fix([index.crpix1,index.crpix2]-0.5) $
                      else sltxy = fix([index.sltpx1ix,index.sltpx2ix]-0.5)
xc = s_psf.xc + (sltxy[0]-s_psf.sltxy[0])
yc = s_psf.yc + (sltxy[1]-s_psf.sltxy[1])

; enforce valid pixel range (just in case, shouldn't be necessary
;                            if missing data was identified properly)
dx1 = dx1 > (min(xc)-3*sz/4)
dx2 = dx2 < (max(xc)+3*sz/4-1)
dy1 = dy1 > (min(yc)-3*sz/4)
dy2 = dy2 < (max(yc)+3*sz/4-1)

; Make sure data range is sufficiently large for width and height of a tile
if dx2-dx1 lt sz-1 or dy2-dy1 lt sz-1 then return,data

; only use PSFs necessary to process image at hand (i.e., when not full frame)
; if no PSFs in range -> return non-deconvolved image
wpsf = where(xc gt dx1 and xc le dx2 and yc gt dy1 and yc le dy2,npsf)
if npsf lt 1 then return,data
psf = s_psf.psf[*,*,wpsf]
xc = xc[wpsf]
yc = yc[wpsf]

; assume grid is regular: find # of grid elements in x and y and their centers
gx = xc[uniq(xc,sort(xc))]
gy = yc[uniq(yc,sort(yc))]
nx = n_elements(gx)
ny = n_elements(gy)
if nx*ny ne npsf then begin
  print,'Error, grid of PSFs not regular'
  stop
endif

; set default for number of RL iterations
if n_elements(niter) ne 1 then niter=0
if niter eq 0 then begin
   d_ave = mean(data[dx1:dx2,dy1:dy2] > 0.0)
   niter = ceil(2.5*d_ave^0.4) > 8       ; not fewer than 8 iterations
endif

; create 4D array for deconvolved tiles
d4d = fltarr(sz,sz,nx,ny)

; loop through grid and process subimages
for k=0,ny-1 do begin
  case k of
  0    : y1 = dy1
  ny-1 : y1 = dy2-sz+1
  else : y1 = gy[k]-sz/2
  endcase
  y2 = y1+sz-1
  for j=0,nx-1 do begin
    wpsf1 = (where(gx[j] eq xc and gy[k] eq yc))[0]
    case j of
    0    : x1 = dx1
    nx-1 : x1 = dx2-sz+1
    else : x1 = gx[j]-sz/2
    endcase
    x2 = x1+sz-1
    psf1 = psf[*,*,wpsf1]
    ; deconvolution with forgiving treatment of missing pixels
    dtil = data[x1:x2,y1:y2]
    ww = where(dtil eq mind,nww)            ; missing data in tile?
    if nww gt maxmis*sz*sz then dtdc=dtil $ ; too much -> use unprocessed
    else begin
       ; preprocess tile if there are missing values
       if nww gt 0 then begin
          dtil[ww] = mean(dtil[where(dtil gt mind)])  ; replace missing w/mean
          dmis = smooth(dtil,9,/edge_trunc) ; smooth out missing data at edges
          dtil[ww] = dmis[ww]               ; smooth out only missing data
       endif
       dtdc = aia_deconvolve_richardsonlucy(dtil,psf1,niter=niter)
       ; postprocess tile if there are bad or missing values
       wnan = where(finite(dtdc,/nan),nnan)
       if nnan gt 0 then dtdc[wnan] = dtil[wnan] ; replace bad output w/input
       if nww  gt 0 then dtdc[ww]  = mind        ; reset missing pixels
    endelse
    d4d[0,0,j,k] = dtdc
  endfor
endfor

; reconstruct tiles into image
out = data < mind ; initialize output array with missing data value
for k=0,ny-1 do begin
  case k of
  0    : begin
           y1  = dy1
	   oy1 = 0
	   oy2 = gy[k]+sz/4-1-y1
	 end
  ny-1 : begin
           y1  = gy[k]-sz/4
	   oy1 = y1-(dy2-sz+1)
	   oy2 = sz-1
	 end
  else : begin
           y1  = gy[k]-sz/4
	   oy1 = sz/4
	   oy2 = 3*sz/4-1
	 end
  endcase
  for j=0,nx-1 do begin
    case j of
    0    : begin
           x1  = dx1
	   ox1 = 0
	   ox2 = gx[j]+sz/4-1-x1
	 end
    nx-1 : begin
           x1  = gx[j]-sz/4
	   ox1 = x1-(dx2-sz+1)
	   ox2 = sz-1
	 end
    else : begin
           x1  = gx[j]-sz/4
	   ox1 = sz/4
	   ox2 = 3*sz/4-1
	 end
    endcase
    out[x1,y1] = d4d[ox1:ox2,oy1:oy2,j,k]
  endfor
endfor

return,out

end