Changeset 473 for trunk/SRC/Interpolation/extrapolate.pro

Timestamp:

02/14/12 16:25:51 (12 years ago)

Author:

smasson

Message:

add options in extrapolation routines

File:

• trunk/SRC/Interpolation/extrapolate.pro (modified) (4 diffs)

trunk/SRC/Interpolation/extrapolate.pro

@@ -64 +64 @@
 ;-
 FUNCTION extrapolate, zinput, maskinput, nb_iteration, FILLXDIR = fillxdir, FILLYDIR = fillydir $
-                      , X_PERIODIC = x_periodic, MINVAL=minval, MAXVAL=maxval, GE0=ge0, _EXTRA=ex
+                      , X_PERIODIC = x_periodic, MINVAL = minval, MAXVAL = maxval, GE0 = ge0 $
+                      , SMWIN = smwin, NSMOOTH = nsmooth, _EXTRA = ex
 ;
   compile_opt idl2, strictarrsubs
@@ -104 +105 @@
     ztmp[temporary(nan)] = 1.e20
     msk = temporary(msk) * temporary(finztmp)
-  ENDIF ELSE finztmp = -1 ; free memory
+  ENDIF ELSE finztmp = -1       ; free memory
   z = temporary(ztmp)
   nx2 = nx+2
@@ -186 +187 @@
       keyword_set(fillydir):zcoast = z[nx2+coast]+z[-nx2+coast]
       ELSE:zcoast = z[1+coast]+z[-1+coast]+z[nx2+coast]+z[-nx2+coast] $
-             +1./sqrt(2)*(z[nx2+1+coast]+z[nx2-1+coast] $
-                          +z[-nx2+1+coast]+z[-nx2-1+coast])
+                    +1./sqrt(2)*(z[nx2+1+coast]+z[nx2-1+coast] $
+                                 +z[-nx2+1+coast]+z[-nx2-1+coast])
     ENDCASE
 ;
@@ -224 +225 @@
 ;
   ENDWHILE
+;
+  z = z[1:nx, 1:ny]
+;
+;---------------------------------------------------------------
+; smooth the filled values
+;---------------------------------------------------------------
+  IF n_elements(smwin) EQ 0 THEN smwin = 0
+  IF smwin EQ 0 OR smwin EQ 1 THEN return, z
+  IF n_elements(nsmooth) EQ 0 THEN nsmooth = 25
+  IF nsmooth EQ 0 THEN return, z
+
+; add extra bands to avoid edge errors...
+
+  one = replicate(1, smwin) 
+  IF keyword_set(x_periodic) THEN BEGIN
+;   nx-smwin   nx-1    0       smwin-1                      nx-smwin      nx-1    0       smwin-1
+;     |----------|     |----------|------------------------------|----------|     |----------|
+;     0      smwin-1  smwin   2*smwin-1                          nx  nx+smwin-1  nx+smwin nx+2*smwin-1 
+    z   = [         z[nx-smwin:nx-1, *],         z,         z[0:smwin-1, *] ]
+    msk = [ maskinput[nx-smwin:nx-1, *], maskinput, maskinput[0:smwin-1, *] ]
+;; IF array_equal(z[0:smwin-1, *], z[nx:nx+smwin-1, *]) NE 1 THEN stop
+;; IF array_equal(z[nx+smwin:nx+2*smwin-1, *], z[smwin:2*smwin-1, *]) NE 1 THEN stop
+  ENDIF ELSE BEGIN 
+    z   = [         one#(z[0, *])[*],         z,         one#(z[nx-1, *])[*] ]
+    msk = [ one#(maskinput[0, *])[*], maskinput, one#(maskinput[nx-1, *])[*] ]
+  ENDELSE 
+  z   = [ [  z[*, 0]#one], [  z], [  z[*, ny-1]#one] ]
+  msk = [ [msk[*, 0]#one], [msk], [msk[*, ny-1]#one] ]
+
+  zorg = z
+  mskm1 = 1-msk
+  FOR i = 1, nsmooth DO BEGIN
+    z = smooth( temporary(z)*mskm1 + zorg*msk, smwin )
+    IF keyword_set(x_periodic) THEN BEGIN
+      z[0:smwin-1, *] = z[nx:nx+smwin-1, *]
+      z[nx+smwin:nx+2*smwin-1, *] =  z[smwin:2*smwin-1, *]
+    ENDIF ELSE BEGIN 
+      z[0:smwin-1, *] = one#(z[smwin, *])[*]
+      z[nx+smwin:nx+2*smwin-1, *] = one#(z[nx+smwin-1, *])[*]
+    ENDELSE 
+    z[*, 0:smwin-1] = z[*, smwin]#one
+    z[*, ny+smwin:ny+2*smwin-1] = z[*, ny+smwin-1]#one
+  ENDFOR  
+  z = temporary(z)*mskm1 + zorg*msk
+  z = (temporary(z))[smwin:nx+smwin-1, smwin:ny+smwin-1]
+;
 ;---------------------------------------------------------------
 ; we return the original size of the array
 ;---------------------------------------------------------------
 ;
-  return, z[1:nx, 1:ny]
+  return, z
 END