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