|
1 #====================================================================== |
|
2 # . N M I N T E R P . L I N E A R |
|
3 # doc: Wed Nov 22 21:01:09 2000 |
|
4 # dlm: Thu Jan 23 13:39:13 2014 |
|
5 # (c) 2000 A.M. Thurnherr |
|
6 # uE-Info: 19 61 NIL 0 0 72 2 2 4 NIL ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # linear interpolation for [fillgaps] |
|
10 |
|
11 # HISTORY: |
|
12 # Apr 6, 2006: - created from [.interp.linear] |
|
13 # Jul 28, 2006: - added xf to ISInit() args |
|
14 # Jul 30, 2006: - added max-gap support |
|
15 # Aug 8, 2008: - documentation |
|
16 # Sep 23, 2011: - added support for %RECNO |
|
17 # - removed xv from interp() args |
|
18 # Oct 19, 2011: - code did not work for %RECNO and trailing missing vals |
|
19 # Jan 23, 2014: - implemented huge speedup for sparse files |
|
20 |
|
21 # NOTES: |
|
22 # - in contrast to the [.interp.*] routines, the [.nminterp.*] routines: |
|
23 # 1) do not assume that x values increase monotonically |
|
24 # 2) can handle nans |
|
25 # 3) cannot be used with multiple buffers |
|
26 # 4) use %RECNO when xfnr==-1 |
|
27 |
|
28 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
29 # |
|
30 # THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading) |
|
31 # |
|
32 # $ISOpts string of allowed options |
|
33 # $ISOptsUsage usage information string for options |
|
34 # |
|
35 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
36 |
|
37 $ISOpts = ""; |
|
38 $ISOptsUsage = ""; |
|
39 |
|
40 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
41 # |
|
42 # &ISUsage() mangle parameters (options, really) |
|
43 # |
|
44 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
45 |
|
46 sub ISUsage() {} |
|
47 |
|
48 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
49 # |
|
50 # &ISInit(f,xf) init interpolation of field f |
|
51 # f field number |
|
52 # xf x field number or -1 for %RECNO |
|
53 # <ret val> none |
|
54 # |
|
55 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
56 |
|
57 sub ISInit($$) {} |
|
58 |
|
59 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
60 # |
|
61 # &interpolate(xf,mg,f,lvr,cr) interpolate field f |
|
62 # xf x field or -1 for %RECNO |
|
63 # mg max allowed x gap |
|
64 # f field number to interpolate |
|
65 # lvr last record with valid y val |
|
66 # cr current record |
|
67 # <ret val> interpolated value |
|
68 # |
|
69 # NB: |
|
70 # - handle f == xf |
|
71 # - return undef if interpolation cannot be carried out |
|
72 # - x VALUES MAY NOT BE MONOTONIC |
|
73 # |
|
74 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
|
75 |
|
76 { my(@nvr); # static: next valid record (per field) |
|
77 |
|
78 sub interpolate($$$$$$) |
|
79 { |
|
80 my($xf,$mg,$f,$lvr,$cr) = @_; |
|
81 |
|
82 return $ants_[$cr][$xf] if ($xf == $f); |
|
83 return $ants_[$lvr][$f] # edge value is ok |
|
84 if ($xf>=0 && equal($ants_[$cr][$xf],$ants_[$lvr][$xf])) || |
|
85 ($xf<0 && $cr==$lvr); |
|
86 |
|
87 unless ($nvr[$f] > $cr) { |
|
88 $nvr[$f] = $cr + 1; |
|
89 while ($nvr[$f] <= $#ants_ && !numberp($ants_[$nvr[$f]][$f])) { $nvr[$f]++; } |
|
90 } |
|
91 |
|
92 return undef if ($nvr[$f] > $#ants_); |
|
93 return undef |
|
94 if ($xf>=0 && abs($ants_[$nvr[$f]][$xf]-$ants_[$lvr][$xf])>$mg); |
|
95 return $ants_[$nvr[$f]][$f] # edge value is ok |
|
96 if ($xf>=0 && equal($ants_[$cr][$xf],$ants_[$nvr[$f]][$xf])) || |
|
97 ($xf<0 && $cr==$nvr[$f]); |
|
98 |
|
99 my($sc) = ($xf >= 0) |
|
100 ? ($ants_[$cr][$xf] - $ants_[$lvr][$xf]) / ($ants_[$nvr[$f]][$xf] - $ants_[$lvr][$xf]) |
|
101 : ($cr - $lvr) / ($nvr[$f] - $lvr); |
|
102 return $ants_[$lvr][$f] + $sc * ($ants_[$nvr[$f]][$f] - $ants_[$lvr][$f]); |
|
103 } |
|
104 } # static scope |
|
105 |
|
106 #---------------------------------------------------------------------- |
|
107 |
|
108 1; |