FUNCTION gridfill, ZIn, nfill=nfill, weightfnc=WeightFnc,status=s ;+ ; NAME: ; gridfill ; PURPOSE: ; Removes empty bins. An empty bin is given the average over all ; neighbours with valid function values. ; CATEGORY: ; Plotting: contours ; CALLING SEQUENCE: ; Z = gridfill(Z,nfile=nFill,weightfnc=WeightFnc,status=S) ; INPUTS: ; Z 2D array; any type ; If of type float, then invalid grid values are marked by !values.f_nan ; OPTIONAL INPUTS: ; nfill=nFill scalar; type integer; default: 0 ; <= 0: fill all empty bins with extrapolated values ; > 0: fill empty bins with more than nFill neighbours ; weightfnc=WeightFnc ; scalar; type string; default: not present ; name of function to calculate weights (see PROCEDURE) ; OUTPUTS: ; Z 2D array; type float ; same as input array with invalid grid values replaced by averages ; (if nfill>0 not all invalid values may have been filled in) ; OPTIONAL OUTPUTS: ; status=S ; 2D array; type byte ; array identifying the extrapolated values: ; = 1 contents of bin is same as valid input value ; = 2 contents of bin is extrapolated value ; nFill>0 only: ; = 0 empty bin with less than nFill neighbours ; Z[.,.] = !values.f_nan (same as input). ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, ArrayLocation ; PROCEDURE: ; > If nFill <= 0 then: ; Step 1: for each empty bin, count the number of non-empty neighbours ; Step 2: find the subset of empty bins with the maximum number of ; non-empty neighbours ; Step 3: for the subset of step 2, calculate the average over the non- ; empty neighbours and assign this average to the empty bin ; Step 4: Go to step 1 ; Repeat until there are no empty bins left. ; ; > nFill > 0: ; Step 1: for each empty bin, count the number of non-empty neighbours ; Step 2: find the subset of empty bins with nFill or more non-empty ; neighbours ; Step 3: for the subset of step 2, calculate the average over the non- ; empty neighbours and assign this average to the empty bin ; Step 4: Return. ; ; > By default neigbouring function values are averaged with weight=1. ; If WeightFnc is specified, it is used to calculate the weights instead. ; The function has the form ; function WeightFnc(I,J,inX,jnY,Z[inX,jnY]) ; I,J are the indices of the bad pixel; inX,jnY are arrays with the ; indices of good neighbour pixels, and Z[inX,jnY] are the ; function values in the good neighbours. ; MODIFICATION HISTORY: ; 1990, Paul Hick (UCSD) ; 24-MAR-1999, Paul Hick (UCSD), converted from Fortran ;- InitVar, nFill, 0 nFill = (nFill > 0) < 8 bNoWeightFnc = n_elements(WeightFnc) EQ 0 Z = float(ZIn) IF (where(finite(Z)))[0] EQ -1 THEN message, 'No valid grid function values' sz = size(Z) nX = sz[1] nY = sz[2] ;------- ; Set up S: PixGood for valid elements of Z; PixBad for invalid elements. ; Then set bad elements in Z to zero PixBad = 0B PixGood = 1B PixAve = 2B S = replicate(PixGood,nX,nY) ; Initialize the status array: all pixels good Bad = where(1B-finite(Z), nBad) IF nBad EQ 0 THEN return, Z ; No empty bins S[Bad] = PixBad ; Set status for bad pixels Side = bytarr(nX,nY,8) Cntr = bytarr(nX,nY) ;------- ; At this point S has value PixGood or PixBad. The bins with value PixBad are to ; be replaced by averaged values. xSide = [-1,0,1,1,1,0,-1,-1] ; Differential indices of neighbours ySide = [1,1,1,0,-1,-1,-1,0] ; Scan through array looking for bad pixels. Identify the bad pixels and ; count the number of valid neighbours Bad = S EQ PixBad ; List of bad pixels (there's at least one) PBad = where(Bad, nPBad) ; Location and # bad pixels ; Intialize the Cntr and Side arrays ; Cntr contains the number of good neigbours for all bad pixels ; Side identifies the good neighbours (1B) and the bad neighbours (0B) ii = nPbad-6 FOR iBad=0L,nPBad-1 DO BEGIN ; Loop over bad pixels P = ArrayLocation(PBad[iBad],size=sz) I = P[0] J = P[1] iX = I+xSide jY = J+ySide ; P identifies neighbours P = where(0 LE iX AND iX LT nX AND 0 LE jY AND jY LT nY) iX = iX[P] jY = jY[P] Good = S[iX,jY] NE PixBad Cntr[I,J ] = round(total(Good)); # good neighbour pixels Side[I,J,P] = Good ; Location of good neighbour pixels ENDFOR REPEAT BEGIN MaxGood = max(Cntr[PBad]) ; Max # of good neighbours (never zero?) ; nFill <= 0: fill in bad pixels with exactly MaxGood good neighbours ; nFill > 0: fill in bad pixels with nFill or more good neighbours ; !! QBad is a subset of PBad QBad = where(Bad AND ( ( nFill EQ 0 AND Cntr EQ MaxGood ) OR $ ( nFill GT 0 AND Cntr GE nFill ) ), nQBad) ; Fill in qualifying pixels with average over good neighbours. FOR iBad=0L,nQBad-1 DO BEGIN P = ArrayLocation(QBad[iBad],size=sz) I = P[0] J = P[1] iX = I+xSide jY = J+ySide P = where(0 LE iX AND iX LT nX AND 0 LE jY AND jY LT nY AND Side[I,J,*]) iX = iX[P] jY = jY[P] CASE bNoWeightFnc OF 1: W = 0.*P+1. 0: W = call_function(Weight,I,J,iX,jY,Z[iX,jY]) ENDCASE Z[I,J] = total(W*Z[iX,jY])/total(W) ; Average of good neighbours ENDFOR ; Update S,nPBad, PBad, Bad, Cntr and Side. ; The loop should not be merged with the previous loop! IF nQBad GT 0 THEN BEGIN ; NOT redundant !! S [QBad] = PixAve ; New status: average value Bad [QBad] = 0B ; QBad pixels not bad anymore Cntr[QBad] = 0B PBad = where(Bad, nPBad); Location and # bad pixels ENDIF FOR iBad=0L,nQBad-1 DO BEGIN ; Loop over pixels just filled with averages P = ArrayLocation(QBad[iBad],size=sz) I = P[0] J = P[1] iX = I+xSide jY = J+ySide P = where(0 LE iX AND iX LT nX AND 0 LE jY AND jY LT nY) iX = iX[P] ; Locations of all neighbours jY = jY[P] ; Look for bad neighbours of the pixel just filled with an average value. ; These bad neighbour pixels now have an additional good neighbour ; (the averaged pixel). RBad = where(Bad[iX,jY],nRBad) IF nRBad GT 0 THEN BEGIN P = P [RBad] iX = iX[RBad] ; Indices of bad neighbours jY = jY[RBad] Cntr[iX,jY] = Cntr[iX,jY]+1B ; One more good neighbour Side[iX,jY,(P+4) mod 8] = 1B ; Opposite neigbhour ENDIF ENDFOR ; Continue filling if nFill=0 ('fill all bad pixels') and there still are ; bad pixels present. ENDREP UNTIL nFill GT 0 OR nPBad EQ 0 RETURN, Z & END