Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3.2/sources/HeuristicLab.ExtLibs/HeuristicLab.ALGLIB/2.3.0/ALGLIB-2.3.0/spline2d.cs

Last change on this file was 2806, checked in by gkronber, 15 years ago

Added plugin for new version of ALGLIB. #875 (Update ALGLIB sources)
File size: 43.6 KB
Line 
1/*************************************************************************
2Copyright (c) 2007, Sergey Bochkanov (ALGLIB project).
3
4>>> SOURCE LICENSE >>>
5This program is free software; you can redistribute it and/or modify
6it under the terms of the GNU General Public License as published by
7the Free Software Foundation (www.fsf.org); either version 2 of the
8License, or (at your option) any later version.
9
10This program is distributed in the hope that it will be useful,
11but WITHOUT ANY WARRANTY; without even the implied warranty of
12MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13GNU General Public License for more details.
14
15A copy of the GNU General Public License is available at
16http://www.fsf.org/licensing/licenses
17
18>>> END OF LICENSE >>>
19*************************************************************************/
20
21using System;
22
23namespace alglib
24{
25    public class spline2d
26    {
27        /*************************************************************************
28        2-dimensional spline inteprolant
29        *************************************************************************/
30        public struct spline2dinterpolant
31        {
32            public int k;
33            public double[] c;
34        };
35
36
37
38
39        public const int spline2dvnum = 12;
40
41
42        /*************************************************************************
43        This subroutine builds bilinear spline coefficients table.
44
45        Input parameters:
46            X   -   spline abscissas, array[0..N-1]
47            Y   -   spline ordinates, array[0..M-1]
48            F   -   function values, array[0..M-1,0..N-1]
49            M,N -   grid size, M>=2, N>=2
50
51        Output parameters:
52            C   -   spline interpolant
53
54          -- ALGLIB PROJECT --
55             Copyright 05.07.2007 by Bochkanov Sergey
56        *************************************************************************/
57        public static void spline2dbuildbilinear(double[] x,
58            double[] y,
59            double[,] f,
60            int m,
61            int n,
62            ref spline2dinterpolant c)
63        {
64            int i = 0;
65            int j = 0;
66            int k = 0;
67            int tblsize = 0;
68            int shift = 0;
69            double t = 0;
70            double[,] dx = new double[0,0];
71            double[,] dy = new double[0,0];
72            double[,] dxy = new double[0,0];
73
74            x = (double[])x.Clone();
75            y = (double[])y.Clone();
76            f = (double[,])f.Clone();
77
78            System.Diagnostics.Debug.Assert(n>=2 & m>=2, "Spline2DBuildBilinear: N<2 or M<2!");
79           
80            //
81            // Sort points
82            //
83            for(j=0; j<=n-1; j++)
84            {
85                k = j;
86                for(i=j+1; i<=n-1; i++)
87                {
88                    if( (double)(x[i])<(double)(x[k]) )
89                    {
90                        k = i;
91                    }
92                }
93                if( k!=j )
94                {
95                    for(i=0; i<=m-1; i++)
96                    {
97                        t = f[i,j];
98                        f[i,j] = f[i,k];
99                        f[i,k] = t;
100                    }
101                    t = x[j];
102                    x[j] = x[k];
103                    x[k] = t;
104                }
105            }
106            for(i=0; i<=m-1; i++)
107            {
108                k = i;
109                for(j=i+1; j<=m-1; j++)
110                {
111                    if( (double)(y[j])<(double)(y[k]) )
112                    {
113                        k = j;
114                    }
115                }
116                if( k!=i )
117                {
118                    for(j=0; j<=n-1; j++)
119                    {
120                        t = f[i,j];
121                        f[i,j] = f[k,j];
122                        f[k,j] = t;
123                    }
124                    t = y[i];
125                    y[i] = y[k];
126                    y[k] = t;
127                }
128            }
129           
130            //
131            // Fill C:
132            //  C[0]            -   length(C)
133            //  C[1]            -   type(C):
134            //                      -1 = bilinear interpolant
135            //                      -3 = general cubic spline
136            //                           (see BuildBicubicSpline)
137            //  C[2]:
138            //      N (x count)
139            //  C[3]:
140            //      M (y count)
141            //  C[4]...C[4+N-1]:
142            //      x[i], i = 0...N-1
143            //  C[4+N]...C[4+N+M-1]:
144            //      y[i], i = 0...M-1
145            //  C[4+N+M]...C[4+N+M+(N*M-1)]:
146            //      f(i,j) table. f(0,0), f(0, 1), f(0,2) and so on...
147            //
148            c.k = 1;
149            tblsize = 4+n+m+n*m;
150            c.c = new double[tblsize-1+1];
151            c.c[0] = tblsize;
152            c.c[1] = -1;
153            c.c[2] = n;
154            c.c[3] = m;
155            for(i=0; i<=n-1; i++)
156            {
157                c.c[4+i] = x[i];
158            }
159            for(i=0; i<=m-1; i++)
160            {
161                c.c[4+n+i] = y[i];
162            }
163            for(i=0; i<=m-1; i++)
164            {
165                for(j=0; j<=n-1; j++)
166                {
167                    shift = i*n+j;
168                    c.c[4+n+m+shift] = f[i,j];
169                }
170            }
171        }
172
173
174        /*************************************************************************
175        This subroutine builds bicubic spline coefficients table.
176
177        Input parameters:
178            X   -   spline abscissas, array[0..N-1]
179            Y   -   spline ordinates, array[0..M-1]
180            F   -   function values, array[0..M-1,0..N-1]
181            M,N -   grid size, M>=2, N>=2
182
183        Output parameters:
184            C   -   spline interpolant
185
186          -- ALGLIB PROJECT --
187             Copyright 05.07.2007 by Bochkanov Sergey
188        *************************************************************************/
189        public static void spline2dbuildbicubic(double[] x,
190            double[] y,
191            double[,] f,
192            int m,
193            int n,
194            ref spline2dinterpolant c)
195        {
196            int i = 0;
197            int j = 0;
198            int k = 0;
199            int tblsize = 0;
200            int shift = 0;
201            double t = 0;
202            double[,] dx = new double[0,0];
203            double[,] dy = new double[0,0];
204            double[,] dxy = new double[0,0];
205
206            x = (double[])x.Clone();
207            y = (double[])y.Clone();
208            f = (double[,])f.Clone();
209
210            System.Diagnostics.Debug.Assert(n>=2 & m>=2, "BuildBicubicSpline: N<2 or M<2!");
211           
212            //
213            // Sort points
214            //
215            for(j=0; j<=n-1; j++)
216            {
217                k = j;
218                for(i=j+1; i<=n-1; i++)
219                {
220                    if( (double)(x[i])<(double)(x[k]) )
221                    {
222                        k = i;
223                    }
224                }
225                if( k!=j )
226                {
227                    for(i=0; i<=m-1; i++)
228                    {
229                        t = f[i,j];
230                        f[i,j] = f[i,k];
231                        f[i,k] = t;
232                    }
233                    t = x[j];
234                    x[j] = x[k];
235                    x[k] = t;
236                }
237            }
238            for(i=0; i<=m-1; i++)
239            {
240                k = i;
241                for(j=i+1; j<=m-1; j++)
242                {
243                    if( (double)(y[j])<(double)(y[k]) )
244                    {
245                        k = j;
246                    }
247                }
248                if( k!=i )
249                {
250                    for(j=0; j<=n-1; j++)
251                    {
252                        t = f[i,j];
253                        f[i,j] = f[k,j];
254                        f[k,j] = t;
255                    }
256                    t = y[i];
257                    y[i] = y[k];
258                    y[k] = t;
259                }
260            }
261           
262            //
263            // Fill C:
264            //  C[0]            -   length(C)
265            //  C[1]            -   type(C):
266            //                      -1 = bilinear interpolant
267            //                           (see BuildBilinearInterpolant)
268            //                      -3 = general cubic spline
269            //  C[2]:
270            //      N (x count)
271            //  C[3]:
272            //      M (y count)
273            //  C[4]...C[4+N-1]:
274            //      x[i], i = 0...N-1
275            //  C[4+N]...C[4+N+M-1]:
276            //      y[i], i = 0...M-1
277            //  C[4+N+M]...C[4+N+M+(N*M-1)]:
278            //      f(i,j) table. f(0,0), f(0, 1), f(0,2) and so on...
279            //  C[4+N+M+N*M]...C[4+N+M+(2*N*M-1)]:
280            //      df(i,j)/dx table.
281            //  C[4+N+M+2*N*M]...C[4+N+M+(3*N*M-1)]:
282            //      df(i,j)/dy table.
283            //  C[4+N+M+3*N*M]...C[4+N+M+(4*N*M-1)]:
284            //      d2f(i,j)/dxdy table.
285            //
286            c.k = 3;
287            tblsize = 4+n+m+4*n*m;
288            c.c = new double[tblsize-1+1];
289            c.c[0] = tblsize;
290            c.c[1] = -3;
291            c.c[2] = n;
292            c.c[3] = m;
293            for(i=0; i<=n-1; i++)
294            {
295                c.c[4+i] = x[i];
296            }
297            for(i=0; i<=m-1; i++)
298            {
299                c.c[4+n+i] = y[i];
300            }
301            bicubiccalcderivatives(ref f, ref x, ref y, m, n, ref dx, ref dy, ref dxy);
302            for(i=0; i<=m-1; i++)
303            {
304                for(j=0; j<=n-1; j++)
305                {
306                    shift = i*n+j;
307                    c.c[4+n+m+shift] = f[i,j];
308                    c.c[4+n+m+n*m+shift] = dx[i,j];
309                    c.c[4+n+m+2*n*m+shift] = dy[i,j];
310                    c.c[4+n+m+3*n*m+shift] = dxy[i,j];
311                }
312            }
313        }
314
315
316        /*************************************************************************
317        This subroutine calculates the value of the bilinear or bicubic spline  at
318        the given point X.
319
320        Input parameters:
321            C   -   coefficients table.
322                    Built by BuildBilinearSpline or BuildBicubicSpline.
323            X, Y-   point
324
325        Result:
326            S(x,y)
327
328          -- ALGLIB PROJECT --
329             Copyright 05.07.2007 by Bochkanov Sergey
330        *************************************************************************/
331        public static double spline2dcalc(ref spline2dinterpolant c,
332            double x,
333            double y)
334        {
335            double result = 0;
336            double v = 0;
337            double vx = 0;
338            double vy = 0;
339            double vxy = 0;
340
341            spline2ddiff(ref c, x, y, ref v, ref vx, ref vy, ref vxy);
342            result = v;
343            return result;
344        }
345
346
347        /*************************************************************************
348        This subroutine calculates the value of the bilinear or bicubic spline  at
349        the given point X and its derivatives.
350
351        Input parameters:
352            C   -   spline interpolant.
353            X, Y-   point
354
355        Output parameters:
356            F   -   S(x,y)
357            FX  -   dS(x,y)/dX
358            FY  -   dS(x,y)/dY
359            FXY -   d2S(x,y)/dXdY
360
361          -- ALGLIB PROJECT --
362             Copyright 05.07.2007 by Bochkanov Sergey
363        *************************************************************************/
364        public static void spline2ddiff(ref spline2dinterpolant c,
365            double x,
366            double y,
367            ref double f,
368            ref double fx,
369            ref double fy,
370            ref double fxy)
371        {
372            int n = 0;
373            int m = 0;
374            double t = 0;
375            double dt = 0;
376            double u = 0;
377            double du = 0;
378            int ix = 0;
379            int iy = 0;
380            int l = 0;
381            int r = 0;
382            int h = 0;
383            int shift1 = 0;
384            int s1 = 0;
385            int s2 = 0;
386            int s3 = 0;
387            int s4 = 0;
388            int sf = 0;
389            int sfx = 0;
390            int sfy = 0;
391            int sfxy = 0;
392            double y1 = 0;
393            double y2 = 0;
394            double y3 = 0;
395            double y4 = 0;
396            double v = 0;
397            double t0 = 0;
398            double t1 = 0;
399            double t2 = 0;
400            double t3 = 0;
401            double u0 = 0;
402            double u1 = 0;
403            double u2 = 0;
404            double u3 = 0;
405
406            System.Diagnostics.Debug.Assert((int)Math.Round(c.c[1])==-1 | (int)Math.Round(c.c[1])==-3, "Spline2DDiff: incorrect C!");
407            n = (int)Math.Round(c.c[2]);
408            m = (int)Math.Round(c.c[3]);
409           
410            //
411            // Binary search in the [ x[0], ..., x[n-2] ] (x[n-1] is not included)
412            //
413            l = 4;
414            r = 4+n-2+1;
415            while( l!=r-1 )
416            {
417                h = (l+r)/2;
418                if( (double)(c.c[h])>=(double)(x) )
419                {
420                    r = h;
421                }
422                else
423                {
424                    l = h;
425                }
426            }
427            t = (x-c.c[l])/(c.c[l+1]-c.c[l]);
428            dt = 1.0/(c.c[l+1]-c.c[l]);
429            ix = l-4;
430           
431            //
432            // Binary search in the [ y[0], ..., y[m-2] ] (y[m-1] is not included)
433            //
434            l = 4+n;
435            r = 4+n+(m-2)+1;
436            while( l!=r-1 )
437            {
438                h = (l+r)/2;
439                if( (double)(c.c[h])>=(double)(y) )
440                {
441                    r = h;
442                }
443                else
444                {
445                    l = h;
446                }
447            }
448            u = (y-c.c[l])/(c.c[l+1]-c.c[l]);
449            du = 1.0/(c.c[l+1]-c.c[l]);
450            iy = l-(4+n);
451           
452            //
453            // Prepare F, dF/dX, dF/dY, d2F/dXdY
454            //
455            f = 0;
456            fx = 0;
457            fy = 0;
458            fxy = 0;
459           
460            //
461            // Bilinear interpolation
462            //
463            if( (int)Math.Round(c.c[1])==-1 )
464            {
465                shift1 = 4+n+m;
466                y1 = c.c[shift1+n*iy+ix];
467                y2 = c.c[shift1+n*iy+(ix+1)];
468                y3 = c.c[shift1+n*(iy+1)+(ix+1)];
469                y4 = c.c[shift1+n*(iy+1)+ix];
470                f = (1-t)*(1-u)*y1+t*(1-u)*y2+t*u*y3+(1-t)*u*y4;
471                fx = (-((1-u)*y1)+(1-u)*y2+u*y3-u*y4)*dt;
472                fy = (-((1-t)*y1)-t*y2+t*y3+(1-t)*y4)*du;
473                fxy = (y1-y2+y3-y4)*du*dt;
474                return;
475            }
476           
477            //
478            // Bicubic interpolation
479            //
480            if( (int)Math.Round(c.c[1])==-3 )
481            {
482               
483                //
484                // Prepare info
485                //
486                t0 = 1;
487                t1 = t;
488                t2 = AP.Math.Sqr(t);
489                t3 = t*t2;
490                u0 = 1;
491                u1 = u;
492                u2 = AP.Math.Sqr(u);
493                u3 = u*u2;
494                sf = 4+n+m;
495                sfx = 4+n+m+n*m;
496                sfy = 4+n+m+2*n*m;
497                sfxy = 4+n+m+3*n*m;
498                s1 = n*iy+ix;
499                s2 = n*iy+(ix+1);
500                s3 = n*(iy+1)+(ix+1);
501                s4 = n*(iy+1)+ix;
502               
503                //
504                // Calculate
505                //
506                v = +(1*c.c[sf+s1]);
507                f = f+v*t0*u0;
508                v = +(1*c.c[sfy+s1]/du);
509                f = f+v*t0*u1;
510                fy = fy+1*v*t0*u0*du;
511                v = -(3*c.c[sf+s1])+3*c.c[sf+s4]-2*c.c[sfy+s1]/du-1*c.c[sfy+s4]/du;
512                f = f+v*t0*u2;
513                fy = fy+2*v*t0*u1*du;
514                v = +(2*c.c[sf+s1])-2*c.c[sf+s4]+1*c.c[sfy+s1]/du+1*c.c[sfy+s4]/du;
515                f = f+v*t0*u3;
516                fy = fy+3*v*t0*u2*du;
517                v = +(1*c.c[sfx+s1]/dt);
518                f = f+v*t1*u0;
519                fx = fx+1*v*t0*u0*dt;
520                v = +(1*c.c[sfxy+s1]/(dt*du));
521                f = f+v*t1*u1;
522                fx = fx+1*v*t0*u1*dt;
523                fy = fy+1*v*t1*u0*du;
524                fxy = fxy+1*v*t0*u0*dt*du;
525                v = -(3*c.c[sfx+s1]/dt)+3*c.c[sfx+s4]/dt-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s4]/(dt*du);
526                f = f+v*t1*u2;
527                fx = fx+1*v*t0*u2*dt;
528                fy = fy+2*v*t1*u1*du;
529                fxy = fxy+2*v*t0*u1*dt*du;
530                v = +(2*c.c[sfx+s1]/dt)-2*c.c[sfx+s4]/dt+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s4]/(dt*du);
531                f = f+v*t1*u3;
532                fx = fx+1*v*t0*u3*dt;
533                fy = fy+3*v*t1*u2*du;
534                fxy = fxy+3*v*t0*u2*dt*du;
535                v = -(3*c.c[sf+s1])+3*c.c[sf+s2]-2*c.c[sfx+s1]/dt-1*c.c[sfx+s2]/dt;
536                f = f+v*t2*u0;
537                fx = fx+2*v*t1*u0*dt;
538                v = -(3*c.c[sfy+s1]/du)+3*c.c[sfy+s2]/du-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s2]/(dt*du);
539                f = f+v*t2*u1;
540                fx = fx+2*v*t1*u1*dt;
541                fy = fy+1*v*t2*u0*du;
542                fxy = fxy+2*v*t1*u0*dt*du;
543                v = +(9*c.c[sf+s1])-9*c.c[sf+s2]+9*c.c[sf+s3]-9*c.c[sf+s4]+6*c.c[sfx+s1]/dt+3*c.c[sfx+s2]/dt-3*c.c[sfx+s3]/dt-6*c.c[sfx+s4]/dt+6*c.c[sfy+s1]/du-6*c.c[sfy+s2]/du-3*c.c[sfy+s3]/du+3*c.c[sfy+s4]/du+4*c.c[sfxy+s1]/(dt*du)+2*c.c[sfxy+s2]/(dt*du)+1*c.c[sfxy+s3]/(dt*du)+2*c.c[sfxy+s4]/(dt*du);
544                f = f+v*t2*u2;
545                fx = fx+2*v*t1*u2*dt;
546                fy = fy+2*v*t2*u1*du;
547                fxy = fxy+4*v*t1*u1*dt*du;
548                v = -(6*c.c[sf+s1])+6*c.c[sf+s2]-6*c.c[sf+s3]+6*c.c[sf+s4]-4*c.c[sfx+s1]/dt-2*c.c[sfx+s2]/dt+2*c.c[sfx+s3]/dt+4*c.c[sfx+s4]/dt-3*c.c[sfy+s1]/du+3*c.c[sfy+s2]/du+3*c.c[sfy+s3]/du-3*c.c[sfy+s4]/du-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s2]/(dt*du)-1*c.c[sfxy+s3]/(dt*du)-2*c.c[sfxy+s4]/(dt*du);
549                f = f+v*t2*u3;
550                fx = fx+2*v*t1*u3*dt;
551                fy = fy+3*v*t2*u2*du;
552                fxy = fxy+6*v*t1*u2*dt*du;
553                v = +(2*c.c[sf+s1])-2*c.c[sf+s2]+1*c.c[sfx+s1]/dt+1*c.c[sfx+s2]/dt;
554                f = f+v*t3*u0;
555                fx = fx+3*v*t2*u0*dt;
556                v = +(2*c.c[sfy+s1]/du)-2*c.c[sfy+s2]/du+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s2]/(dt*du);
557                f = f+v*t3*u1;
558                fx = fx+3*v*t2*u1*dt;
559                fy = fy+1*v*t3*u0*du;
560                fxy = fxy+3*v*t2*u0*dt*du;
561                v = -(6*c.c[sf+s1])+6*c.c[sf+s2]-6*c.c[sf+s3]+6*c.c[sf+s4]-3*c.c[sfx+s1]/dt-3*c.c[sfx+s2]/dt+3*c.c[sfx+s3]/dt+3*c.c[sfx+s4]/dt-4*c.c[sfy+s1]/du+4*c.c[sfy+s2]/du+2*c.c[sfy+s3]/du-2*c.c[sfy+s4]/du-2*c.c[sfxy+s1]/(dt*du)-2*c.c[sfxy+s2]/(dt*du)-1*c.c[sfxy+s3]/(dt*du)-1*c.c[sfxy+s4]/(dt*du);
562                f = f+v*t3*u2;
563                fx = fx+3*v*t2*u2*dt;
564                fy = fy+2*v*t3*u1*du;
565                fxy = fxy+6*v*t2*u1*dt*du;
566                v = +(4*c.c[sf+s1])-4*c.c[sf+s2]+4*c.c[sf+s3]-4*c.c[sf+s4]+2*c.c[sfx+s1]/dt+2*c.c[sfx+s2]/dt-2*c.c[sfx+s3]/dt-2*c.c[sfx+s4]/dt+2*c.c[sfy+s1]/du-2*c.c[sfy+s2]/du-2*c.c[sfy+s3]/du+2*c.c[sfy+s4]/du+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s2]/(dt*du)+1*c.c[sfxy+s3]/(dt*du)+1*c.c[sfxy+s4]/(dt*du);
567                f = f+v*t3*u3;
568                fx = fx+3*v*t2*u3*dt;
569                fy = fy+3*v*t3*u2*du;
570                fxy = fxy+9*v*t2*u2*dt*du;
571                return;
572            }
573        }
574
575
576        /*************************************************************************
577        This subroutine unpacks two-dimensional spline into the coefficients table
578
579        Input parameters:
580            C   -   spline interpolant.
581
582        Result:
583            M, N-   grid size (x-axis and y-axis)
584            Tbl -   coefficients table, unpacked format,
585                    [0..(N-1)*(M-1)-1, 0..19].
586                    For I = 0...M-2, J=0..N-2:
587                        K =  I*(N-1)+J
588                        Tbl[K,0] = X[j]
589                        Tbl[K,1] = X[j+1]
590                        Tbl[K,2] = Y[i]
591                        Tbl[K,3] = Y[i+1]
592                        Tbl[K,4] = C00
593                        Tbl[K,5] = C01
594                        Tbl[K,6] = C02
595                        Tbl[K,7] = C03
596                        Tbl[K,8] = C10
597                        Tbl[K,9] = C11
598                        ...
599                        Tbl[K,19] = C33
600                    On each grid square spline is equals to:
601                        S(x) = SUM(c[i,j]*(x^i)*(y^j), i=0..3, j=0..3)
602                        t = x-x[j]
603                        u = y-y[i]
604
605          -- ALGLIB PROJECT --
606             Copyright 29.06.2007 by Bochkanov Sergey
607        *************************************************************************/
608        public static void spline2dunpack(ref spline2dinterpolant c,
609            ref int m,
610            ref int n,
611            ref double[,] tbl)
612        {
613            int i = 0;
614            int j = 0;
615            int ci = 0;
616            int cj = 0;
617            int k = 0;
618            int p = 0;
619            int shift = 0;
620            int s1 = 0;
621            int s2 = 0;
622            int s3 = 0;
623            int s4 = 0;
624            int sf = 0;
625            int sfx = 0;
626            int sfy = 0;
627            int sfxy = 0;
628            double y1 = 0;
629            double y2 = 0;
630            double y3 = 0;
631            double y4 = 0;
632            double dt = 0;
633            double du = 0;
634
635            System.Diagnostics.Debug.Assert((int)Math.Round(c.c[1])==-3 | (int)Math.Round(c.c[1])==-1, "SplineUnpack2D: incorrect C!");
636            n = (int)Math.Round(c.c[2]);
637            m = (int)Math.Round(c.c[3]);
638            tbl = new double[(n-1)*(m-1)-1+1, 19+1];
639           
640            //
641            // Fill
642            //
643            for(i=0; i<=m-2; i++)
644            {
645                for(j=0; j<=n-2; j++)
646                {
647                    p = i*(n-1)+j;
648                    tbl[p,0] = c.c[4+j];
649                    tbl[p,1] = c.c[4+j+1];
650                    tbl[p,2] = c.c[4+n+i];
651                    tbl[p,3] = c.c[4+n+i+1];
652                    dt = 1/(tbl[p,1]-tbl[p,0]);
653                    du = 1/(tbl[p,3]-tbl[p,2]);
654                   
655                    //
656                    // Bilinear interpolation
657                    //
658                    if( (int)Math.Round(c.c[1])==-1 )
659                    {
660                        for(k=4; k<=19; k++)
661                        {
662                            tbl[p,k] = 0;
663                        }
664                        shift = 4+n+m;
665                        y1 = c.c[shift+n*i+j];
666                        y2 = c.c[shift+n*i+(j+1)];
667                        y3 = c.c[shift+n*(i+1)+(j+1)];
668                        y4 = c.c[shift+n*(i+1)+j];
669                        tbl[p,4] = y1;
670                        tbl[p,4+1*4+0] = y2-y1;
671                        tbl[p,4+0*4+1] = y4-y1;
672                        tbl[p,4+1*4+1] = y3-y2-y4+y1;
673                    }
674                   
675                    //
676                    // Bicubic interpolation
677                    //
678                    if( (int)Math.Round(c.c[1])==-3 )
679                    {
680                        sf = 4+n+m;
681                        sfx = 4+n+m+n*m;
682                        sfy = 4+n+m+2*n*m;
683                        sfxy = 4+n+m+3*n*m;
684                        s1 = n*i+j;
685                        s2 = n*i+(j+1);
686                        s3 = n*(i+1)+(j+1);
687                        s4 = n*(i+1)+j;
688                        tbl[p,4+0*4+0] = +(1*c.c[sf+s1]);
689                        tbl[p,4+0*4+1] = +(1*c.c[sfy+s1]/du);
690                        tbl[p,4+0*4+2] = -(3*c.c[sf+s1])+3*c.c[sf+s4]-2*c.c[sfy+s1]/du-1*c.c[sfy+s4]/du;
691                        tbl[p,4+0*4+3] = +(2*c.c[sf+s1])-2*c.c[sf+s4]+1*c.c[sfy+s1]/du+1*c.c[sfy+s4]/du;
692                        tbl[p,4+1*4+0] = +(1*c.c[sfx+s1]/dt);
693                        tbl[p,4+1*4+1] = +(1*c.c[sfxy+s1]/(dt*du));
694                        tbl[p,4+1*4+2] = -(3*c.c[sfx+s1]/dt)+3*c.c[sfx+s4]/dt-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s4]/(dt*du);
695                        tbl[p,4+1*4+3] = +(2*c.c[sfx+s1]/dt)-2*c.c[sfx+s4]/dt+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s4]/(dt*du);
696                        tbl[p,4+2*4+0] = -(3*c.c[sf+s1])+3*c.c[sf+s2]-2*c.c[sfx+s1]/dt-1*c.c[sfx+s2]/dt;
697                        tbl[p,4+2*4+1] = -(3*c.c[sfy+s1]/du)+3*c.c[sfy+s2]/du-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s2]/(dt*du);
698                        tbl[p,4+2*4+2] = +(9*c.c[sf+s1])-9*c.c[sf+s2]+9*c.c[sf+s3]-9*c.c[sf+s4]+6*c.c[sfx+s1]/dt+3*c.c[sfx+s2]/dt-3*c.c[sfx+s3]/dt-6*c.c[sfx+s4]/dt+6*c.c[sfy+s1]/du-6*c.c[sfy+s2]/du-3*c.c[sfy+s3]/du+3*c.c[sfy+s4]/du+4*c.c[sfxy+s1]/(dt*du)+2*c.c[sfxy+s2]/(dt*du)+1*c.c[sfxy+s3]/(dt*du)+2*c.c[sfxy+s4]/(dt*du);
699                        tbl[p,4+2*4+3] = -(6*c.c[sf+s1])+6*c.c[sf+s2]-6*c.c[sf+s3]+6*c.c[sf+s4]-4*c.c[sfx+s1]/dt-2*c.c[sfx+s2]/dt+2*c.c[sfx+s3]/dt+4*c.c[sfx+s4]/dt-3*c.c[sfy+s1]/du+3*c.c[sfy+s2]/du+3*c.c[sfy+s3]/du-3*c.c[sfy+s4]/du-2*c.c[sfxy+s1]/(dt*du)-1*c.c[sfxy+s2]/(dt*du)-1*c.c[sfxy+s3]/(dt*du)-2*c.c[sfxy+s4]/(dt*du);
700                        tbl[p,4+3*4+0] = +(2*c.c[sf+s1])-2*c.c[sf+s2]+1*c.c[sfx+s1]/dt+1*c.c[sfx+s2]/dt;
701                        tbl[p,4+3*4+1] = +(2*c.c[sfy+s1]/du)-2*c.c[sfy+s2]/du+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s2]/(dt*du);
702                        tbl[p,4+3*4+2] = -(6*c.c[sf+s1])+6*c.c[sf+s2]-6*c.c[sf+s3]+6*c.c[sf+s4]-3*c.c[sfx+s1]/dt-3*c.c[sfx+s2]/dt+3*c.c[sfx+s3]/dt+3*c.c[sfx+s4]/dt-4*c.c[sfy+s1]/du+4*c.c[sfy+s2]/du+2*c.c[sfy+s3]/du-2*c.c[sfy+s4]/du-2*c.c[sfxy+s1]/(dt*du)-2*c.c[sfxy+s2]/(dt*du)-1*c.c[sfxy+s3]/(dt*du)-1*c.c[sfxy+s4]/(dt*du);
703                        tbl[p,4+3*4+3] = +(4*c.c[sf+s1])-4*c.c[sf+s2]+4*c.c[sf+s3]-4*c.c[sf+s4]+2*c.c[sfx+s1]/dt+2*c.c[sfx+s2]/dt-2*c.c[sfx+s3]/dt-2*c.c[sfx+s4]/dt+2*c.c[sfy+s1]/du-2*c.c[sfy+s2]/du-2*c.c[sfy+s3]/du+2*c.c[sfy+s4]/du+1*c.c[sfxy+s1]/(dt*du)+1*c.c[sfxy+s2]/(dt*du)+1*c.c[sfxy+s3]/(dt*du)+1*c.c[sfxy+s4]/(dt*du);
704                    }
705                   
706                    //
707                    // Rescale Cij
708                    //
709                    for(ci=0; ci<=3; ci++)
710                    {
711                        for(cj=0; cj<=3; cj++)
712                        {
713                            tbl[p,4+ci*4+cj] = tbl[p,4+ci*4+cj]*Math.Pow(dt, ci)*Math.Pow(du, cj);
714                        }
715                    }
716                }
717            }
718        }
719
720
721        /*************************************************************************
722        This subroutine performs linear transformation of the spline argument.
723
724        Input parameters:
725            C       -   spline interpolant
726            AX, BX  -   transformation coefficients: x = A*t + B
727            AY, BY  -   transformation coefficients: y = A*u + B
728        Result:
729            C   -   transformed spline
730
731          -- ALGLIB PROJECT --
732             Copyright 30.06.2007 by Bochkanov Sergey
733        *************************************************************************/
734        public static void spline2dlintransxy(ref spline2dinterpolant c,
735            double ax,
736            double bx,
737            double ay,
738            double by)
739        {
740            int i = 0;
741            int j = 0;
742            int n = 0;
743            int m = 0;
744            double v = 0;
745            double[] x = new double[0];
746            double[] y = new double[0];
747            double[,] f = new double[0,0];
748            int typec = 0;
749
750            typec = (int)Math.Round(c.c[1]);
751            System.Diagnostics.Debug.Assert(typec==-3 | typec==-1, "Spline2DLinTransXY: incorrect C!");
752            n = (int)Math.Round(c.c[2]);
753            m = (int)Math.Round(c.c[3]);
754            x = new double[n-1+1];
755            y = new double[m-1+1];
756            f = new double[m-1+1, n-1+1];
757            for(j=0; j<=n-1; j++)
758            {
759                x[j] = c.c[4+j];
760            }
761            for(i=0; i<=m-1; i++)
762            {
763                y[i] = c.c[4+n+i];
764            }
765            for(i=0; i<=m-1; i++)
766            {
767                for(j=0; j<=n-1; j++)
768                {
769                    f[i,j] = c.c[4+n+m+i*n+j];
770                }
771            }
772           
773            //
774            // Special case: AX=0 or AY=0
775            //
776            if( (double)(ax)==(double)(0) )
777            {
778                for(i=0; i<=m-1; i++)
779                {
780                    v = spline2dcalc(ref c, bx, y[i]);
781                    for(j=0; j<=n-1; j++)
782                    {
783                        f[i,j] = v;
784                    }
785                }
786                if( typec==-3 )
787                {
788                    spline2dbuildbicubic(x, y, f, m, n, ref c);
789                }
790                if( typec==-1 )
791                {
792                    spline2dbuildbilinear(x, y, f, m, n, ref c);
793                }
794                ax = 1;
795                bx = 0;
796            }
797            if( (double)(ay)==(double)(0) )
798            {
799                for(j=0; j<=n-1; j++)
800                {
801                    v = spline2dcalc(ref c, x[j], by);
802                    for(i=0; i<=m-1; i++)
803                    {
804                        f[i,j] = v;
805                    }
806                }
807                if( typec==-3 )
808                {
809                    spline2dbuildbicubic(x, y, f, m, n, ref c);
810                }
811                if( typec==-1 )
812                {
813                    spline2dbuildbilinear(x, y, f, m, n, ref c);
814                }
815                ay = 1;
816                by = 0;
817            }
818           
819            //
820            // General case: AX<>0, AY<>0
821            // Unpack, scale and pack again.
822            //
823            for(j=0; j<=n-1; j++)
824            {
825                x[j] = (x[j]-bx)/ax;
826            }
827            for(i=0; i<=m-1; i++)
828            {
829                y[i] = (y[i]-by)/ay;
830            }
831            if( typec==-3 )
832            {
833                spline2dbuildbicubic(x, y, f, m, n, ref c);
834            }
835            if( typec==-1 )
836            {
837                spline2dbuildbilinear(x, y, f, m, n, ref c);
838            }
839        }
840
841
842        /*************************************************************************
843        This subroutine performs linear transformation of the spline.
844
845        Input parameters:
846            C   -   spline interpolant.
847            A, B-   transformation coefficients: S2(x,y) = A*S(x,y) + B
848           
849        Output parameters:
850            C   -   transformed spline
851
852          -- ALGLIB PROJECT --
853             Copyright 30.06.2007 by Bochkanov Sergey
854        *************************************************************************/
855        public static void spline2dlintransf(ref spline2dinterpolant c,
856            double a,
857            double b)
858        {
859            int i = 0;
860            int j = 0;
861            int n = 0;
862            int m = 0;
863            double[] x = new double[0];
864            double[] y = new double[0];
865            double[,] f = new double[0,0];
866            int typec = 0;
867
868            typec = (int)Math.Round(c.c[1]);
869            System.Diagnostics.Debug.Assert(typec==-3 | typec==-1, "Spline2DLinTransXY: incorrect C!");
870            n = (int)Math.Round(c.c[2]);
871            m = (int)Math.Round(c.c[3]);
872            x = new double[n-1+1];
873            y = new double[m-1+1];
874            f = new double[m-1+1, n-1+1];
875            for(j=0; j<=n-1; j++)
876            {
877                x[j] = c.c[4+j];
878            }
879            for(i=0; i<=m-1; i++)
880            {
881                y[i] = c.c[4+n+i];
882            }
883            for(i=0; i<=m-1; i++)
884            {
885                for(j=0; j<=n-1; j++)
886                {
887                    f[i,j] = a*c.c[4+n+m+i*n+j]+b;
888                }
889            }
890            if( typec==-3 )
891            {
892                spline2dbuildbicubic(x, y, f, m, n, ref c);
893            }
894            if( typec==-1 )
895            {
896                spline2dbuildbilinear(x, y, f, m, n, ref c);
897            }
898        }
899
900
901        /*************************************************************************
902        This subroutine makes the copy of the spline model.
903
904        Input parameters:
905            C   -   spline interpolant
906
907        Output parameters:
908            CC  -   spline copy
909
910          -- ALGLIB PROJECT --
911             Copyright 29.06.2007 by Bochkanov Sergey
912        *************************************************************************/
913        public static void spline2dcopy(ref spline2dinterpolant c,
914            ref spline2dinterpolant cc)
915        {
916            int n = 0;
917            int i_ = 0;
918
919            System.Diagnostics.Debug.Assert(c.k==1 | c.k==3, "Spline2DCopy: incorrect C!");
920            cc.k = c.k;
921            n = (int)Math.Round(c.c[0]);
922            cc.c = new double[n];
923            for(i_=0; i_<=n-1;i_++)
924            {
925                cc.c[i_] = c.c[i_];
926            }
927        }
928
929
930        /*************************************************************************
931        Serialization of the spline interpolant
932
933        INPUT PARAMETERS:
934            B   -   spline interpolant
935
936        OUTPUT PARAMETERS:
937            RA      -   array of real numbers which contains interpolant,
938                        array[0..RLen-1]
939            RLen    -   RA lenght
940
941          -- ALGLIB --
942             Copyright 17.08.2009 by Bochkanov Sergey
943        *************************************************************************/
944        public static void spline2dserialize(ref spline2dinterpolant c,
945            ref double[] ra,
946            ref int ralen)
947        {
948            int clen = 0;
949            int i_ = 0;
950            int i1_ = 0;
951
952            System.Diagnostics.Debug.Assert(c.k==1 | c.k==3, "Spline2DSerialize: incorrect C!");
953            clen = (int)Math.Round(c.c[0]);
954            ralen = 3+clen;
955            ra = new double[ralen];
956            ra[0] = ralen;
957            ra[1] = spline2dvnum;
958            ra[2] = c.k;
959            i1_ = (0) - (3);
960            for(i_=3; i_<=3+clen-1;i_++)
961            {
962                ra[i_] = c.c[i_+i1_];
963            }
964        }
965
966
967        /*************************************************************************
968        Unserialization of the spline interpolant
969
970        INPUT PARAMETERS:
971            RA  -   array of real numbers which contains interpolant,
972
973        OUTPUT PARAMETERS:
974            B   -   spline interpolant
975
976          -- ALGLIB --
977             Copyright 17.08.2009 by Bochkanov Sergey
978        *************************************************************************/
979        public static void spline2dunserialize(ref double[] ra,
980            ref spline2dinterpolant c)
981        {
982            int clen = 0;
983            int i_ = 0;
984            int i1_ = 0;
985
986            System.Diagnostics.Debug.Assert((int)Math.Round(ra[1])==spline2dvnum, "Spline2DUnserialize: corrupted array!");
987            c.k = (int)Math.Round(ra[2]);
988            clen = (int)Math.Round(ra[3]);
989            c.c = new double[clen];
990            i1_ = (3) - (0);
991            for(i_=0; i_<=clen-1;i_++)
992            {
993                c.c[i_] = ra[i_+i1_];
994            }
995        }
996
997
998        /*************************************************************************
999        Bicubic spline resampling
1000
1001        Input parameters:
1002            A           -   function values at the old grid,
1003                            array[0..OldHeight-1, 0..OldWidth-1]
1004            OldHeight   -   old grid height, OldHeight>1
1005            OldWidth    -   old grid width, OldWidth>1
1006            NewHeight   -   new grid height, NewHeight>1
1007            NewWidth    -   new grid width, NewWidth>1
1008           
1009        Output parameters:
1010            B           -   function values at the new grid,
1011                            array[0..NewHeight-1, 0..NewWidth-1]
1012
1013          -- ALGLIB routine --
1014             15 May, 2007
1015             Copyright by Bochkanov Sergey
1016        *************************************************************************/
1017        public static void spline2dresamplebicubic(ref double[,] a,
1018            int oldheight,
1019            int oldwidth,
1020            ref double[,] b,
1021            int newheight,
1022            int newwidth)
1023        {
1024            double[,] buf = new double[0,0];
1025            double[] x = new double[0];
1026            double[] y = new double[0];
1027            spline1d.spline1dinterpolant c = new spline1d.spline1dinterpolant();
1028            int i = 0;
1029            int j = 0;
1030            int mw = 0;
1031            int mh = 0;
1032
1033            System.Diagnostics.Debug.Assert(oldwidth>1 & oldheight>1, "Spline2DResampleBicubic: width/height less than 1");
1034            System.Diagnostics.Debug.Assert(newwidth>1 & newheight>1, "Spline2DResampleBicubic: width/height less than 1");
1035           
1036            //
1037            // Prepare
1038            //
1039            mw = Math.Max(oldwidth, newwidth);
1040            mh = Math.Max(oldheight, newheight);
1041            b = new double[newheight-1+1, newwidth-1+1];
1042            buf = new double[oldheight-1+1, newwidth-1+1];
1043            x = new double[Math.Max(mw, mh)-1+1];
1044            y = new double[Math.Max(mw, mh)-1+1];
1045           
1046            //
1047            // Horizontal interpolation
1048            //
1049            for(i=0; i<=oldheight-1; i++)
1050            {
1051               
1052                //
1053                // Fill X, Y
1054                //
1055                for(j=0; j<=oldwidth-1; j++)
1056                {
1057                    x[j] = (double)(j)/((double)(oldwidth-1));
1058                    y[j] = a[i,j];
1059                }
1060               
1061                //
1062                // Interpolate and place result into temporary matrix
1063                //
1064                spline1d.spline1dbuildcubic(x, y, oldwidth, 0, 0.0, 0, 0.0, ref c);
1065                for(j=0; j<=newwidth-1; j++)
1066                {
1067                    buf[i,j] = spline1d.spline1dcalc(ref c, (double)(j)/((double)(newwidth-1)));
1068                }
1069            }
1070           
1071            //
1072            // Vertical interpolation
1073            //
1074            for(j=0; j<=newwidth-1; j++)
1075            {
1076               
1077                //
1078                // Fill X, Y
1079                //
1080                for(i=0; i<=oldheight-1; i++)
1081                {
1082                    x[i] = (double)(i)/((double)(oldheight-1));
1083                    y[i] = buf[i,j];
1084                }
1085               
1086                //
1087                // Interpolate and place result into B
1088                //
1089                spline1d.spline1dbuildcubic(x, y, oldheight, 0, 0.0, 0, 0.0, ref c);
1090                for(i=0; i<=newheight-1; i++)
1091                {
1092                    b[i,j] = spline1d.spline1dcalc(ref c, (double)(i)/((double)(newheight-1)));
1093                }
1094            }
1095        }
1096
1097
1098        /*************************************************************************
1099        Bilinear spline resampling
1100
1101        Input parameters:
1102            A           -   function values at the old grid,
1103                            array[0..OldHeight-1, 0..OldWidth-1]
1104            OldHeight   -   old grid height, OldHeight>1
1105            OldWidth    -   old grid width, OldWidth>1
1106            NewHeight   -   new grid height, NewHeight>1
1107            NewWidth    -   new grid width, NewWidth>1
1108
1109        Output parameters:
1110            B           -   function values at the new grid,
1111                            array[0..NewHeight-1, 0..NewWidth-1]
1112
1113          -- ALGLIB routine --
1114             09.07.2007
1115             Copyright by Bochkanov Sergey
1116        *************************************************************************/
1117        public static void spline2dresamplebilinear(ref double[,] a,
1118            int oldheight,
1119            int oldwidth,
1120            ref double[,] b,
1121            int newheight,
1122            int newwidth)
1123        {
1124            int i = 0;
1125            int j = 0;
1126            int l = 0;
1127            int c = 0;
1128            double t = 0;
1129            double u = 0;
1130
1131            b = new double[newheight-1+1, newwidth-1+1];
1132            for(i=0; i<=newheight-1; i++)
1133            {
1134                for(j=0; j<=newwidth-1; j++)
1135                {
1136                    l = i*(oldheight-1)/(newheight-1);
1137                    if( l==oldheight-1 )
1138                    {
1139                        l = oldheight-2;
1140                    }
1141                    u = (double)(i)/((double)(newheight-1))*(oldheight-1)-l;
1142                    c = j*(oldwidth-1)/(newwidth-1);
1143                    if( c==oldwidth-1 )
1144                    {
1145                        c = oldwidth-2;
1146                    }
1147                    t = (double)(j*(oldwidth-1))/((double)(newwidth-1))-c;
1148                    b[i,j] = (1-t)*(1-u)*a[l,c]+t*(1-u)*a[l,c+1]+t*u*a[l+1,c+1]+(1-t)*u*a[l+1,c];
1149                }
1150            }
1151        }
1152
1153
1154        /*************************************************************************
1155        Internal subroutine.
1156        Calculation of the first derivatives and the cross-derivative.
1157        *************************************************************************/
1158        private static void bicubiccalcderivatives(ref double[,] a,
1159            ref double[] x,
1160            ref double[] y,
1161            int m,
1162            int n,
1163            ref double[,] dx,
1164            ref double[,] dy,
1165            ref double[,] dxy)
1166        {
1167            int i = 0;
1168            int j = 0;
1169            double[] xt = new double[0];
1170            double[] ft = new double[0];
1171            double[] c = new double[0];
1172            double s = 0;
1173            double ds = 0;
1174            double d2s = 0;
1175
1176            dx = new double[m-1+1, n-1+1];
1177            dy = new double[m-1+1, n-1+1];
1178            dxy = new double[m-1+1, n-1+1];
1179           
1180            //
1181            // dF/dX
1182            //
1183            xt = new double[n-1+1];
1184            ft = new double[n-1+1];
1185            for(i=0; i<=m-1; i++)
1186            {
1187                for(j=0; j<=n-1; j++)
1188                {
1189                    xt[j] = x[j];
1190                    ft[j] = a[i,j];
1191                }
1192                spline3.buildcubicspline(xt, ft, n, 0, 0.0, 0, 0.0, ref c);
1193                for(j=0; j<=n-1; j++)
1194                {
1195                    spline3.splinedifferentiation(ref c, x[j], ref s, ref ds, ref d2s);
1196                    dx[i,j] = ds;
1197                }
1198            }
1199           
1200            //
1201            // dF/dY
1202            //
1203            xt = new double[m-1+1];
1204            ft = new double[m-1+1];
1205            for(j=0; j<=n-1; j++)
1206            {
1207                for(i=0; i<=m-1; i++)
1208                {
1209                    xt[i] = y[i];
1210                    ft[i] = a[i,j];
1211                }
1212                spline3.buildcubicspline(xt, ft, m, 0, 0.0, 0, 0.0, ref c);
1213                for(i=0; i<=m-1; i++)
1214                {
1215                    spline3.splinedifferentiation(ref c, y[i], ref s, ref ds, ref d2s);
1216                    dy[i,j] = ds;
1217                }
1218            }
1219           
1220            //
1221            // d2F/dXdY
1222            //
1223            xt = new double[n-1+1];
1224            ft = new double[n-1+1];
1225            for(i=0; i<=m-1; i++)
1226            {
1227                for(j=0; j<=n-1; j++)
1228                {
1229                    xt[j] = x[j];
1230                    ft[j] = dy[i,j];
1231                }
1232                spline3.buildcubicspline(xt, ft, n, 0, 0.0, 0, 0.0, ref c);
1233                for(j=0; j<=n-1; j++)
1234                {
1235                    spline3.splinedifferentiation(ref c, x[j], ref s, ref ds, ref d2s);
1236                    dxy[i,j] = ds;
1237                }
1238            }
1239        }
1240    }
1241}
Note: See TracBrowser for help on using the repository browser.