des.math.method.approx.interp source code

1 module des.math.method.approx.interp;
2 
3 import std.algorithm;
4 import std.exception;
5 import std.math;
6 import std.traits;
7 
8 version(unittest) import des.math.linear.vector;
9 
10 import des.ts;
11 import des.math.combin;
12 import des.stdx.traits : hasBasicMathOp;
13 
14 ///
15 struct InterpolateTableData(T) if( hasBasicMathOp!T ) { float key; T val; }
16 
17 ///
18 auto lineInterpolate(T)( in InterpolateTableData!T[] tbl, float k, bool line_end=false )
19     if( hasBasicMathOp!T )
20 {
21     enforce( tbl.length > 1 );
22     size_t i = tbl.length - find!"a.key > b"( tbl, k ).length;
23     if( !line_end )
24     {
25         if( i <= 0 ) return tbl[0].val;
26         else if( i >= tbl.length ) return tbl[$-1].val;
27     }
28     else
29     {
30         if( i < 1 ) i = 1;
31         else if( i > tbl.length-1 ) i = tbl.length-1;
32     }
33 
34     auto a = tbl[i-1];
35     auto b = tbl[i];
36     return cast(T)( a.val + ( b.val - a.val ) * ( ( k - a.key ) / ( b.key - a.key ) ) );
37 }
38 
39 ///
40 unittest
41 {
42     alias InterpolateTableData!float TT;
43     auto tbl =
44         [
45         TT( 0, 10 ),
46         TT( 10, 18 ),
47         TT( 25, 20 ),
48         TT( 50, 13 ),
49         TT( 55, 25 )
50         ];
51 
52     assertEq( lineInterpolate( tbl, 0 ), 10 );
53     assertEq( lineInterpolate( tbl, 5 ), 14 );
54     assertEq( lineInterpolate( tbl, 10 ), 18 );
55     assertEq( lineInterpolate( tbl, -10 ), 10 );
56     assertEq( lineInterpolate( tbl, 80 ), 25 );
57 }
58 
59 unittest
60 {
61     alias InterpolateTableData!double TD;
62     auto tbl =
63         [
64         TD( 0, 0 ),
65         TD( 1, 1 ),
66         TD( 2, 3 ),
67         TD( 3, 4 )
68         ];
69     assertEq( lineInterpolate( tbl, 5, true ), 6 );
70     assertEq( lineInterpolate( tbl, -3, true ), -3 );
71 }
72 
73 unittest
74 {
75     alias InterpolateTableData!vec3 TC;
76     auto tbl =
77         [
78         TC( 0, vec3(1,0,0) ),
79         TC( 1, vec3(0,1,0) ),
80         TC( 2, vec3(0,0,1) )
81         ];
82 
83     assertEq( lineInterpolate( tbl, -1 ), vec3(1,0,0) );
84     assertEq( lineInterpolate( tbl, 0 ), vec3(1,0,0) );
85     assertEq( lineInterpolate( tbl, 0.5 ), vec3(0.5,0.5,0) );
86     assertEq( lineInterpolate( tbl, 3 ), vec3(0,0,1) );
87 }
88 
89 @property bool canBezierInterpolate(T,F)()
90 { return is( typeof( T.init * F.init + T.init * F.init ) == T ) && isNumeric!F; }
91 
92 pure nothrow auto bezierInterpolation(T,F=float)( in T[] pts, F t )
93 if( canBezierInterpolate!(T,F) )
94 in
95 {
96     assert( t >= 0 && t <= 1 );
97     assert( pts.length > 0 );
98 }
99 body
100 {
101     auto N = pts.length-1;
102     auto omt = 1.0 - t;
103     T res = pts[0] * pow(omt,N) + pts[$-1] * pow(t,N);
104     foreach( i; 1 .. N )
105         res = res + pts[i] * pow(t,i) * pow(omt,N-i) * combinationCount(N,i);
106     return res;
107 }
108 
109 unittest
110 {
111     auto pts = [ vec2(0,0), vec2(2,2), vec2(4,0) ];
112     assertEq( bezierInterpolation( pts, 0.5 ), vec2(2,1) );
113 }
114 
115 pure nothrow auto fixBezierSpline(T,F=float)( in T[] pts, size_t steps )
116 {
117     auto step = cast(F)(1.0) / cast(F)(steps-1);
118     auto res = new T[](steps);
119     for( auto i=0; i < steps; i++ )
120         res[i] = bezierInterpolation( pts, step * i );
121     return res;
122 }
123 
124 interface BezierSplineInterpolator(T,F=float)
125 if( canBezierInterpolate!(T,F) )
126 { T[] opCall( in T[] ); }
127 
128 class FixStepsBezierSplineInterpolator(T,F=float) : BezierSplineInterpolator!(T,F)
129 {
130     size_t steps;
131     this( size_t steps ) { this.steps = steps; }
132 
133     T[] opCall( in T[] pts )
134     { return fixBezierSpline!(T,F)( pts, steps ); }
135 }
136 
137 unittest
138 {
139     enum size_t len = 100;
140     auto fbi = new FixStepsBezierSplineInterpolator!(vec2)(len);
141     auto pts = [ vec2(0,0), vec2(2,2), vec2(4,0) ];
142     auto res = fbi( pts );
143     assertEq( res.length, len );
144 }
145 
146 /+ функция критеря должна быть функцией хевисада +/
147 auto criteriaBezierSpline(T,F=float)( in T[] pts, bool delegate(in T[], in T) criteria, F min_step=1e-5 )
148 if( canBezierInterpolate!(T,F) )
149 in
150 {
151     assert( pts.length > 1 );
152     assert( criteria !is null );
153     assert( min_step > 0 );
154     assert( min_step < cast(F)(.25) / cast(F)(pts.length-1) );
155 }
156 body
157 {
158     F step = cast(F)(.25) / cast(F)(pts.length-1);
159 
160     auto ret = [ bezierInterpolation( pts, 0 ), bezierInterpolation( pts, min_step ) ];
161 
162     auto t = min_step;
163     auto o = step;
164 
165     while( true )
166     {
167         if( 1-t < o ) o = 1-t;
168         auto buf = bezierInterpolation( pts, t+o );
169 
170         if( criteria( ret, buf ) )
171         {
172             t += o;
173             if( t >= 1 ) return ret ~ bezierInterpolation( pts, 1 );
174             continue;
175         }
176         else
177         {
178             o /= 2.0;
179             if( o > min_step ) continue;
180             else t += o;
181         }
182 
183         buf = bezierInterpolation( pts, t );
184 
185         if( t > 1 ) return ret ~ bezierInterpolation( pts, 1 );
186 
187         ret ~= buf;
188         o = step;
189     }
190 }
191 
192 /+ функция критеря должна быть функцией хевисада +/
193 auto filterBezierSpline(T,F=float)( in T[] pts, bool delegate(in T[], in T) criteria )
194 if( canBezierInterpolate!(T,F) )
195 in
196 {
197     assert( pts.length > 1 );
198     assert( criteria !is null );
199 }
200 body
201 {
202     auto ret = [ pts[0] ];
203 
204     for( auto i = 1; i < pts.length; i++ )
205     {
206         while( i < pts.length && criteria(ret, pts[i]) ) i++;
207 
208         if( ret[$-1] == pts[i-1] )
209             ret ~= pts[i];
210         else
211             ret ~= pts[i-1];
212     }
213 
214     if( ret[$-1] != pts[$-1] )
215         ret ~= pts[$-1];
216 
217     return ret;
218 }
219 
220 auto angleSplineCriteria(T)( float angle )
221     if( isSpecVector!(2,float,T) || isSpecVector!(3,float,T) )
222 {
223     auto cos_angle = cos( angle );
224     return ( in T[] accepted, in T newpoint )
225     {
226         if( accepted.length < 2 ) return false;
227 
228         auto cc = [ accepted[$-2], accepted[$-1], newpoint ];
229 
230         auto a = (cc[1]-cc[0]).e;
231         auto b = (cc[2]-cc[1]).e;
232 
233         return dot(a,b) >= cos_angle;
234     };
235 }
236 
237 auto lengthSplineCriteria(T)( float len )
238 if( isSpecVector!(2,float,T) || isSpecVector!(3,float,T) )
239 in{ assert( len > 0 ); } body
240 {
241     return ( in T[] accepted, in T newpoint )
242     { return (accepted[$-1] - newpoint).len2 <= len*len; };
243 }
244 
245 unittest
246 {
247     auto pts = [ vec2(0,0), vec2(5,2), vec2(-1,2), vec2(4,0) ];
248     auto pp = criteriaBezierSpline( pts, angleSplineCriteria!vec2(0.05) );
249 
250     assert( pp.length > pts.length );
251     assertEq( pp[0], pts[0] );
252     assertEq( pp[$-1], pts[$-1] );
253 }
254 
255 unittest
256 {
257     auto pts = [ vec2(0,0), vec2(5,2), vec2(-1,2), vec2(4,0) ];
258     auto pp = filterBezierSpline( fixBezierSpline( pts, 1000 ), lengthSplineCriteria!vec2(0.05) );
259 
260     assert( pp.length > pts.length );
261     assertEq( pp[0], pts[0] );
262     assertEq( pp[$-1], pts[$-1] );
263 }