NOAA-PMEL / Ferret

The Ferret program from NOAA/PMEL
https://ferret.pmel.noaa.gov/Ferret/
The Unlicense
55 stars 20 forks source link

Smoothing transforms that don't peel back the edges #626

Open karlmsmith opened 6 years ago

karlmsmith commented 6 years ago

Reported by steven.c.hankin on 2 Sep 2005 21:31 UTC Ferret's smoothing transforms require a "window" of data. As a result, they are unable to return an answer within 1/2 window of the edge of a domain. How about if parallel to

@SBX box smoothed
@SBN binomial smoothed @SWL Welch smoothed @SHN Hanning smoothed @SPZ Parzen smoothed

we also offered all of the same names, but with an "e" (edges) attached:

@SBXe box smoothed
@SBNe binomial smoothed @SWLe Welch smoothed @SHNe Hanning smoothed @SPZe Parzen smoothed

that included the right mathematics to "optimillay" fill/smooth the edge points, too?

Migrated-From: http://dunkel.pmel.noaa.gov/trac/ferret/ticket/1337

karlmsmith commented 6 years ago

Comment by @AnsleyManke on 6 Sep 2005 16:24 UTC Here are the suggestions from Ryo Furue furue`@hawaii`.edu on how we might do the filling; email of 9/2/2005

Hello Ferret users,

I noticed that smoothing transformations (@SPZ and friends) leave endpoints undefined. For example,

yes? let pi = 4 atan(1) yes? let func sin(pix[x=-1:1:0.2]) yes? plot/line/symbol/hlimits=-1.2:1.2 func yes? plot/line/symbol/hlimits=-1.2:1.2/ov func[x=@spz]

Even though func has datapoints at the edges (x = -1 and x = 1), func[x=@spz] has missing values there. I'm wondering if those transformations can be "improved", at least in simple cases.

Let me take @SPZ as an example. @SPZ is a 1-2-1 weighted moving average:

g(i) = [f(i-1) + 2 f(i) + f(i+1)] / 4

where f is defined for i = 1, 2, . . ., N. What to do with g(1) and g(N), for which f(0) and f(N+1) would be required?

One solution is to define g(1) == f(1) and g(N) == f(N). But, I think a better solution is to make the smoothing behave like diffusion with a no-flux boundary condition(*). That is, we define extra gridpoints at i = 0 and i = N + 1 and define the values of f(i) there as f(0) == f(1) and f(N+1) == f(N). And then we compute g(1) and g(N) using the formula above:

g(1) = [f(0) + 2 f(1) + f(2)] / 4
     = [f(1) + 2 f(1) + f(2)] / 4
     = [       3 f(1) + f(2)] / 4

and likewise for g(N).

We could apply the same no-flux condition around "interior" missing datapoints, too.

Advantages of this approach are that averages are conserved: g[i=1:N@sum] == f[i=1:N@sum] (This property corresponds to the conservation of heat in the diffusion equation with no-flux boundary condition.), and that it matches our intuitive understanding of "mixing" neighboring values.

I haven't given a serious thought to other lowpass filters; the above consideration may or may not apply to them.

A practical reason for proposing this is that I don't like losing vectors or shading near the boundaries when smoothing is applied.

Regards, Ryo

(*)In fact, the 1-2-1 moving average can be cast into the diffusion equation:

g(i) = [f(i-1) + 2 f(i) + f(i+1)] / 4 <--> g(i) - f(i) = [f(i-1) - 2 f(i) + f(i+1)] / 4 <--> [g(i) - f(i)] / delta_t = kappa [f(i-1) - 2 f(i) + f(i+1)] / (delta_x)2

with a suitable choice of delta_t, kappa, and delta_x. This is a finite-difference form of the diffusion equation with g(i) being the value of the next timestep.

I think the 1-4-6-4-1 smoothing can be likewise cast into a diffusion equation with biharmonic diffusion.