des.math.method.approx.interp source code

1 /+
2 The MIT License (MIT)
3 
4     Copyright (c) <2013> <Oleg Butko (deviator), Anton Akzhigitov (Akzwar)>
5 
6     Permission is hereby granted, free of charge, to any person obtaining a copy
7     of this software and associated documentation files (the "Software"), to deal
8     in the Software without restriction, including without limitation the rights
9     to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10     copies of the Software, and to permit persons to whom the Software is
11     furnished to do so, subject to the following conditions:
12 
13     The above copyright notice and this permission notice shall be included in
14     all copies or substantial portions of the Software.
15 
16     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17     IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18     FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19     AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20     LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21     OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
22     THE SOFTWARE.
23 +/
24 
25 module des.math.method.approx.interp;
26 
27 import std.algorithm;
28 import std.exception;
29 import std.math;
30 import std.traits;
31 
32 version(unittest) import des.math.linear.vector;
33 
34 import des.math.combin;
35 public import des.math.basic.traits;
36 
37 ///
38 struct InterpolateTableData(T) if( hasBasicMathOp!T ) { float key; T val; }
39 
40 ///
41 auto lineInterpolate(T)( in InterpolateTableData!T[] tbl, float k, bool line_end=false )
42     if( hasBasicMathOp!T )
43 {
44     enforce( tbl.length > 1 );
45     size_t i = tbl.length - find!"a.key > b"( tbl, k ).length;
46     if( !line_end )
47     {
48         if( i <= 0 ) return tbl[0].val;
49         else if( i >= tbl.length ) return tbl[$-1].val;
50     }
51     else
52     {
53         if( i < 1 ) i = 1;
54         else if( i > tbl.length-1 ) i = tbl.length-1;
55     }
56 
57     auto a = tbl[i-1];
58     auto b = tbl[i];
59     return cast(T)( a.val + ( b.val - a.val ) * ( ( k - a.key ) / ( b.key - a.key ) ) );
60 }
61 
62 ///
63 unittest
64 {
65     alias InterpolateTableData!float TT;
66     auto tbl =
67         [
68         TT( 0, 10 ),
69         TT( 10, 18 ),
70         TT( 25, 20 ),
71         TT( 50, 13 ),
72         TT( 55, 25 )
73         ];
74 
75     assert( lineInterpolate( tbl, 0 ) == 10 );
76     assert( lineInterpolate( tbl, 5 ) == 14 );
77     assert( lineInterpolate( tbl, 10 ) == 18 );
78     assert( lineInterpolate( tbl, -10 ) == 10 );
79     assert( lineInterpolate( tbl, 80 ) == 25 );
80 }
81 
82 unittest
83 {
84     alias InterpolateTableData!double TD;
85     auto tbl =
86         [
87         TD( 0, 0 ),
88         TD( 1, 1 ),
89         TD( 2, 3 ),
90         TD( 3, 4 )
91         ];
92     assert( lineInterpolate( tbl, 5, true ) == 6 );
93     assert( lineInterpolate( tbl, -3, true ) == -3 );
94 }
95 
96 unittest
97 {
98     alias InterpolateTableData!col3 TC;
99     auto tbl =
100         [
101         TC( 0, col3(1,0,0) ),
102         TC( 1, col3(0,1,0) ),
103         TC( 2, col3(0,0,1) )
104         ];
105 
106     assert( lineInterpolate( tbl, -1 ) == col3(1,0,0) );
107     assert( lineInterpolate( tbl, 0 ) == col3(1,0,0) );
108     assert( lineInterpolate( tbl, 0.5 ) == col3(0.5,0.5,0) );
109     assert( lineInterpolate( tbl, 3 ) == col3(0,0,1) );
110 }
111 
112 @property bool canBezierInterpolate(T,F)()
113 { return is( typeof( T.init * F.init + T.init * F.init ) == T ) && isNumeric!F; }
114 
115 pure nothrow auto bezierInterpolation(T,F=float)( in T[] pts, F t )
116 if( canBezierInterpolate!(T,F) )
117 in
118 {
119     assert( t >= 0 && t <= 1 );
120     assert( pts.length > 0 );
121 }
122 body
123 {
124     auto N = pts.length-1;
125     auto omt = 1.0 - t;
126     T res = pts[0] * pow(omt,N) + pts[$-1] * pow(t,N);
127     for( auto i=1; i < N; i++ )
128         res = res + pts[i] * pow(t,i) * pow(omt,N-i) * combination(N,i);
129     return res;
130 }
131 
132 unittest
133 {
134     auto pts = [ vec2(0,0), vec2(2,2), vec2(4,0) ];
135     assert( bezierInterpolation( pts, 0.5 ) == vec2(2,1) );
136 }
137 
138 pure nothrow auto fixBezierSpline(T,F=float)( in T[] pts, size_t steps )
139 {
140     auto step = cast(F)(1.0) / cast(F)(steps-1);
141     auto res = new T[](steps);
142     for( auto i=0; i < steps; i++ )
143         res[i] = bezierInterpolation( pts, step * i );
144     return res;
145 }
146 
147 interface BezierSplineInterpolator(T,F=float)
148 if( canBezierInterpolate!(T,F) )
149 { T[] opCall( in T[] ); }
150 
151 class FixStepsBezierSplineInterpolator(T,F=float) : BezierSplineInterpolator!(T,F)
152 {
153     size_t steps;
154     this( size_t steps ) { this.steps = steps; }
155 
156     T[] opCall( in T[] pts )
157     { return fixBezierSpline!(T,F)( pts, steps ); }
158 }
159 
160 unittest
161 {
162     enum size_t len = 100;
163     auto fbi = new FixStepsBezierSplineInterpolator!(vec2)(len);
164     auto pts = [ vec2(0,0), vec2(2,2), vec2(4,0) ];
165     auto res = fbi( pts );
166     assert( res.length == len );
167 }
168 
169 /+ функция критеря должна быть функцией хевисада +/
170 auto criteriaBezierSpline(T,F=float)( in T[] pts, bool delegate(in T[], in T) criteria, F min_step=1e-5 )
171 if( canBezierInterpolate!(T,F) )
172 in
173 {
174     assert( pts.length > 1 );
175     assert( criteria !is null );
176     assert( min_step > 0 );
177     assert( min_step < cast(F)(.25) / cast(F)(pts.length-1) );
178 }
179 body
180 {
181     F step = cast(F)(.25) / cast(F)(pts.length-1);
182 
183     auto ret = [ bezierInterpolation( pts, 0 ), bezierInterpolation( pts, min_step ) ];
184 
185     auto t = min_step;
186     auto o = step;
187 
188     while( true )
189     {
190         if( 1-t < o ) o = 1-t;
191         auto buf = bezierInterpolation( pts, t+o );
192 
193         if( criteria( ret, buf ) )
194         {
195             t += o;
196             if( t >= 1 ) return ret ~ bezierInterpolation( pts, 1 );
197             continue;
198         }
199         else
200         {
201             o /= 2.0;
202             if( o > min_step ) continue;
203             else t += o;
204         }
205 
206         buf = bezierInterpolation( pts, t );
207 
208         if( t > 1 ) return ret ~ bezierInterpolation( pts, 1 );
209 
210         ret ~= buf;
211         o = step;
212     }
213 }
214 
215 /+ функция критеря должна быть функцией хевисада +/
216 auto filterBezierSpline(T,F=float)( in T[] pts, bool delegate(in T[], in T) criteria )
217 if( canBezierInterpolate!(T,F) )
218 in
219 {
220     assert( pts.length > 1 );
221     assert( criteria !is null );
222 }
223 body
224 {
225     auto ret = [ pts[0] ];
226 
227     for( auto i = 1; i < pts.length; i++ )
228     {
229         while( i < pts.length && criteria(ret, pts[i]) ) i++;
230 
231         if( ret[$-1] == pts[i-1] )
232             ret ~= pts[i];
233         else
234             ret ~= pts[i-1];
235     }
236 
237     if( ret[$-1] != pts[$-1] )
238         ret ~= pts[$-1];
239 
240     return ret;
241 }
242 
243 auto angleSplineCriteria(T)( float angle )
244     if( isCompatibleVector!(2,float,T) || isCompatibleVector!(3,float,T) )
245 {
246     auto cos_angle = cos( angle );
247     return ( in T[] accepted, in T newpoint )
248     {
249         if( accepted.length < 2 ) return false;
250 
251         auto cc = [ accepted[$-2], accepted[$-1], newpoint ];
252 
253         auto a = (cc[1]-cc[0]).e;
254         auto b = (cc[2]-cc[1]).e;
255 
256         return dot(a,b) >= cos_angle;
257     };
258 }
259 
260 auto lengthSplineCriteria(T)( float len )
261 if( isCompatibleVector!(2,float,T) || isCompatibleVector!(3,float,T) )
262 in{ assert( len > 0 ); } body
263 {
264     return ( in T[] accepted, in T newpoint )
265     { return (accepted[$-1] - newpoint).len2 <= len*len; };
266 }
267 
268 unittest
269 {
270     auto pts = [ vec2(0,0), vec2(5,2), vec2(-1,2), vec2(4,0) ];
271     auto pp = criteriaBezierSpline( pts, angleSplineCriteria!vec2(0.05) );
272 
273     assert( pp.length > pts.length );
274     assert( pp[0] == pts[0] );
275     assert( pp[$-1] == pts[$-1] );
276 }
277 
278 unittest
279 {
280     auto pts = [ vec2(0,0), vec2(5,2), vec2(-1,2), vec2(4,0) ];
281     auto pp = filterBezierSpline( fixBezierSpline( pts, 1000 ), lengthSplineCriteria!vec2(0.05) );
282 
283     assert( pp.length > pts.length );
284     assert( pp[0] == pts[0] );
285     assert( pp[$-1] == pts[$-1] );
286 }