Module: ClimbFactor::Filtering

Defined in:
lib/climb_factor/filtering.rb

Class Method Summary collapse

Class Method Details

.do(v0, w) ⇒ Object 



37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
 
# File 'lib/climb_factor/filtering.rb', line 37

def self.do(v0,w)
  # inputs:
  #   v0 is a list of floating-point numbers, representing the value of a function at equally spaced points on the x axis
  #   w = width of rectangular window to convolve with; will be made even if it isn't
  # output:
  #   a fresh array in the same format as v0, doesn't modify v0
  # v0's length should be a power of 2 and >=2

  if w<=1 then return v0 end

  if w%2==1 then w=w+1 end

  # remove DC and detrend, so that start and end are both at 0
  #        -- https://www.dsprelated.com/showthread/comp.dsp/175408-1.php
  # After filtering, we put these back in.
  # v0 = original, which we don't touch
  # v = detrended
  # v1 = filtered
  # For the current method, it's only necessary that n is even, not a power of 2, and detrending
  # isn't actually needed.
  v = v0.dup
  n = v.length
  slope = (v.last-v[0])/(n.to_f-1.0)
  c = -v[0]
  0.upto(n-1) { |i|
    v[i] = v[i] - (c + slope*i)
  }

  # Copy the unfiltered data over as a default. On the initial and final portions, where part of the
  # rectangular kernel hangs over the end, we don't attempt to do any filtering. Using the filter
  # on those portions, even with appropriate normalization, would bias the (x,y) points, effectively
  # moving the start and finish line inward.
  v1 = v.dup

  # convolve with a rectangle of width w:
  sum = 0
  count = 0
  # Sum the initial portion for use in the average for the first filtered data point:
  sum_left = 0.0
  0.upto(w-1) { |i|
    break if i>n/2-1 # this happens in the unusual case where w isn't less than n; we're guaranteed that n is even
    sum_left = sum_left+v[i]
  }
  # the filter is applied to the middle portion, from w to n-w:
  if w<n then
    sum = sum_left
    w.upto(n) { |i|
      if i>=v.length then break end # wasn't part of original algorithm, but needed for some data sets...?
      if v[i].nil?
        raise("coordinate #{i} of #{w}..#{n} is nil, length of vector is #{v.length}, for v0=#{v0.length}, w=#{w}")
      end
      sum = sum + v[i]-v[i-w]
      j = i-w/2
      break if j>n-w
      if j>=w && j<=n-w then
        v1[j] = sum/w
      end
    }
  end

  # To avoid a huge discontinuity in the elevation when the filter turns on, turn it on gradually
  # in the initial and final segments of length w:
  # FIXME: leaves a small discontinuity
  sum_left = 0.0
  sum_right = 0.0
  nn = 0
  0.upto(2*w+1) { |i|
    break if i>n/2-1 # unusual case, see above
    j = n-i-1
    sum_left = sum_left+v[i]
    sum_right = sum_right+v[j]
    nn = nn+1
    if i%2==0 then
      ii = i/2
      jj = n-i/2-1
      v1[ii] = sum_left/nn
      v1[jj] = sum_right/nn
    end
  }

  # put DC and trend back in:
  0.upto(n-1) { |i|
    v1[i] = v1[i] + (c + slope*i)
  }
  return v1
end

.resample_and_filter_hv(hv, filtering, target_resampling_resolution: 30.0) ⇒ Object 



5
6
7
8
9
10
11
12
13
14
15
16
17
18
 
# File 'lib/climb_factor/filtering.rb', line 5

def self.resample_and_filter_hv(hv,filtering,target_resampling_resolution:30.0)
  # The default for target_resampling_resolution is chosen because typically public DEM data isn't any better than this horizontal resolution anyway.
  if filtering>target_resampling_resolution then filtering=target_resampling_resolution end
  k = (filtering/target_resampling_resolution).floor # Set resampling resolution so that it divides filtering.
  while target_resampling_resolution*k<filtering do k+=1 end
  if k%2==1 then k+=1 end # make it even
  resampling_resolution = filtering/k
  resampled_v = resample_hv(hv,resampling_resolution)
  result = []
  0.upto(resampled_v.length-1) { |i|
    result.push([resampling_resolution*i,resampled_v[i]])
  }
  return result
end

.resample_hv(hv, desired_h_resolution) ⇒ Object 



20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 
# File 'lib/climb_factor/filtering.rb', line 20

def self.resample_hv(hv,desired_h_resolution)
  # Inputs a list of [h,v] points, returns a list of v values interpolated to achieve the given constant horizontal resolution (or slightly more).
  # This is similar to add_resolution_and_check_size_limit(), but that routine kludgingly does multiple things.
  tot_h = hv.last[0]-hv[0][0]
  n = (tot_h/desired_h_resolution).ceil+1
  dh = tot_h/(n-1) # actual resolution, which will be a tiny bit better
  result = []
  m = 0 # index into hv
  0.upto(n-1) { |i| # i is an index into resampled output
    h = hv[0][0]+dh*i
    while m<hv.length-1 && hv[m+1][0]<h do m += 1 end
    s = (h-hv[m][0])/(hv[m+1][0]-hv[m][0]) # interpolation fraction ranging from 0 to 1
    result.push(CfMath.linear_interp(hv[m][1],hv[m+1][1],s))
  }
  return result
end