1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Application Programming Interface //
9 // //
10 // Library: SAGA_API //
11 // //
12 //-------------------------------------------------------//
13 // //
14 // mat_spline.cpp //
15 // //
16 // Copyright (C) 2005 by Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. //
22 // //
23 // This library is free software; you can redistribute //
24 // it and/or modify it under the terms of the GNU Lesser //
25 // General Public License as published by the Free //
26 // Software Foundation, either version 2.1 of the //
27 // License, or (at your option) any later version. //
28 // //
29 // This library is distributed in the hope that it will //
30 // be useful, but WITHOUT ANY WARRANTY; without even the //
31 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
32 // PARTICULAR PURPOSE. See the GNU Lesser General Public //
33 // License for more details. //
34 // //
35 // You should have received a copy of the GNU Lesser //
36 // General Public License along with this program; if //
37 // not, see <http://www.gnu.org/licenses/>. //
38 // //
39 //-------------------------------------------------------//
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Goettingen //
44 // Goldschmidtstr. 5 //
45 // 37077 Goettingen //
46 // Germany //
47 // //
48 // e-mail: oconrad@saga-gis.org //
49 // //
50 ///////////////////////////////////////////////////////////
51
52 //---------------------------------------------------------
53 #include "mat_tools.h"
54
55
56 ///////////////////////////////////////////////////////////
57 // //
58 // //
59 // //
60 ///////////////////////////////////////////////////////////
61
62 //---------------------------------------------------------
CSG_Spline(void)63 CSG_Spline::CSG_Spline(void)
64 {
65 m_bCreated = false;
66 }
67
68 //---------------------------------------------------------
~CSG_Spline(void)69 CSG_Spline::~CSG_Spline(void)
70 {
71 Destroy();
72 }
73
74
75 ///////////////////////////////////////////////////////////
76 // //
77 ///////////////////////////////////////////////////////////
78
79 //---------------------------------------------------------
Destroy(void)80 void CSG_Spline::Destroy(void)
81 {
82 m_x.Destroy();
83 m_y.Destroy();
84 m_z.Destroy();
85
86 m_bCreated = false;
87 }
88
89 //---------------------------------------------------------
Create(double * xValues,double * yValues,int nValues,double yA,double yB)90 bool CSG_Spline::Create(double *xValues, double *yValues, int nValues, double yA, double yB)
91 {
92 Destroy();
93
94 for(int i=0; i<nValues; i++)
95 {
96 Add(xValues[i], yValues[i]);
97 }
98
99 return( _Create(yA, yB) );
100 }
101
102 //---------------------------------------------------------
Create(double yA,double yB)103 bool CSG_Spline::Create(double yA, double yB)
104 {
105 return( _Create(yA, yB) );
106 }
107
108
109 ///////////////////////////////////////////////////////////
110 // //
111 ///////////////////////////////////////////////////////////
112
113 //---------------------------------------------------------
Add(double x,double y)114 void CSG_Spline::Add(double x, double y)
115 {
116 m_bCreated = false;
117
118 m_x.Add_Row(x);
119 m_y.Add_Row(y);
120 }
121
122
123 ///////////////////////////////////////////////////////////
124 // //
125 ///////////////////////////////////////////////////////////
126
127 //---------------------------------------------------------
_Create(double yA,double yB)128 bool CSG_Spline::_Create(double yA, double yB)
129 {
130 if( Get_Count() < 3 )
131 {
132 return( false );
133 }
134
135 //-----------------------------------------------------
136 int i, k, n = Get_Count();
137 double p, qn, sig, un;
138 CSG_Vector u;
139
140 //-----------------------------------------------------
141 CSG_Index Index(n, m_x.Get_Data());
142 CSG_Vector x(m_x), y(m_y);
143
144 for(i=0; i<n; i++)
145 {
146 m_x[i] = x[Index[i]];
147 m_y[i] = y[Index[i]];
148 }
149
150 //-----------------------------------------------------
151 u .Create(n);
152 m_z.Create(n);
153
154 if( yA > 0.99e30 )
155 {
156 m_z[0] = u[0] = 0.;
157 }
158 else
159 {
160 m_z[0] = -0.5;
161 u [0] = (3. / (m_x[1] - m_x[0])) * ((m_y[1] - m_y[0]) / (m_x[1] - m_x[0]) - yA);
162 }
163
164 //-----------------------------------------------------
165 for(i=1; i<n-1; i++)
166 {
167 sig = (m_x[i] - m_x[i - 1]) / (m_x[i + 1] - m_x[i - 1]);
168 p = sig * m_z[i - 1] + 2.;
169 m_z[i] = (sig - 1.) / p;
170 u [i] = (m_y[i + 1] - m_y[i ]) / (m_x[i + 1] - m_x[i ])
171 - (m_y[i ] - m_y[i - 1]) / (m_x[i ] - m_x[i - 1]);
172 u [i] = (6. * u[i] / (m_x[i + 1] - m_x[i - 1]) - sig * u[i - 1]) / p;
173 }
174
175 if( yB > 0.99e30 )
176 {
177 qn = un = 0.;
178 }
179 else
180 {
181 qn = 0.5;
182 un = (3. / (m_x[n - 1] - m_x[n - 2]))
183 * (yB - (m_y[n - 1] - m_y[n - 2])
184 / (m_x[n - 1] - m_x[n - 2]));
185 }
186
187 m_z[n - 1] = (un - qn * u[n - 2]) / (qn * m_z[n - 2] + 1.);
188
189 for(k=n-2; k>=0; k--)
190 {
191 m_z[k] = m_z[k] * m_z[k + 1] + u[k];
192 }
193
194 //-----------------------------------------------------
195 m_bCreated = true;
196
197 return( true );
198 }
199
200
201 ///////////////////////////////////////////////////////////
202 // //
203 ///////////////////////////////////////////////////////////
204
205 //---------------------------------------------------------
Get_Value(double x,double & y)206 bool CSG_Spline::Get_Value(double x, double &y)
207 {
208 if( m_bCreated || Create() )
209 {
210 int klo, khi, k;
211 double h, b, a;
212
213 klo = 0;
214 khi = Get_Count() - 1;
215
216 while( khi - klo > 1 )
217 {
218 k = (khi+klo) >> 1;
219
220 if( m_x[k] > x )
221 {
222 khi = k;
223 }
224 else
225 {
226 klo = k;
227 }
228 }
229
230 h = m_x[khi] - m_x[klo];
231
232 if( h != 0. )
233 {
234 a = (m_x[khi] - x) / h;
235 b = (x - m_x[klo]) / h;
236
237 y = a * m_y[klo] + b * m_y[khi]
238 + ((a*a*a - a) * m_z[klo] + (b*b*b - b) * m_z[khi]) * (h*h) / 6.;
239
240 return( true );
241 }
242 }
243
244 return( false );
245 }
246
247 //---------------------------------------------------------
Get_Value(double x)248 double CSG_Spline::Get_Value(double x)
249 {
250 Get_Value(x, x);
251
252 return( x );
253 }
254
255
256 ///////////////////////////////////////////////////////////
257 // //
258 // //
259 // //
260 ///////////////////////////////////////////////////////////
261
262 //---------------------------------------------------------
263 //
264 // Based on:
265 //
266 // Thin Plate Spline demo/example in C++
267 // Copyright (C) 2003, 2005 by Jarno Elonen
268 //
269 // Permission to use, copy, modify, distribute and sell this software
270 // and its documentation for any purpose is hereby granted without fee,
271 // provided that the above copyright notice appear in all copies and
272 // that both that copyright notice and this permission notice appear
273 // in supporting documentation. The authors make no representations
274 // about the suitability of this software for any purpose.
275 // It is provided "as is" without express or implied warranty.
276 //
277 // Reference:
278 // - Donato, G., Belongie, S. (2002):
279 // 'Approximation Methods for Thin Plate Spline Mappings and Principal Warps'
280 //
281 //---------------------------------------------------------
282
283 ///////////////////////////////////////////////////////////
284 // //
285 // //
286 // //
287 ///////////////////////////////////////////////////////////
288
289 //---------------------------------------------------------
CSG_Thin_Plate_Spline(void)290 CSG_Thin_Plate_Spline::CSG_Thin_Plate_Spline(void)
291 {
292 }
293
294 //---------------------------------------------------------
~CSG_Thin_Plate_Spline(void)295 CSG_Thin_Plate_Spline::~CSG_Thin_Plate_Spline(void)
296 {
297 Destroy();
298 }
299
300
301 ///////////////////////////////////////////////////////////
302 // //
303 ///////////////////////////////////////////////////////////
304
305 //---------------------------------------------------------
Destroy(void)306 bool CSG_Thin_Plate_Spline::Destroy(void)
307 {
308 m_Points.Clear();
309 m_V.Destroy();
310
311 return( true );
312 }
313
314
315 ///////////////////////////////////////////////////////////
316 // //
317 ///////////////////////////////////////////////////////////
318
319 //---------------------------------------------------------
_Get_hDistance(TSG_Point_Z A,TSG_Point_Z B)320 double CSG_Thin_Plate_Spline::_Get_hDistance(TSG_Point_Z A, TSG_Point_Z B)
321 {
322 A.x -= B.x;
323 A.y -= B.y;
324
325 return( sqrt(A.x*A.x + A.y*A.y) );
326 }
327
328 //---------------------------------------------------------
_Get_Base_Funtion(double x)329 double CSG_Thin_Plate_Spline::_Get_Base_Funtion(double x)
330 {
331 return( x > 0. ? x*x * log(x) : 0. );
332 }
333
334 //---------------------------------------------------------
_Get_Base_Funtion(TSG_Point_Z A,double x,double y)335 double CSG_Thin_Plate_Spline::_Get_Base_Funtion(TSG_Point_Z A, double x, double y)
336 {
337 x -= A.x;
338 y -= A.y;
339
340 double d = sqrt(x*x + y*y);
341
342 return( d > 0. ? d*d * log(d) : 0. );
343 }
344
345
346 ///////////////////////////////////////////////////////////
347 // //
348 ///////////////////////////////////////////////////////////
349
350 //---------------------------------------------------------
351 // Calculate Thin Plate Spline (TPS) weights from control points.
352 //
Create(double Regularization,bool bSilent)353 bool CSG_Thin_Plate_Spline::Create(double Regularization, bool bSilent)
354 {
355 bool bResult = false;
356 int n;
357 CSG_Matrix M;
358
359 //-----------------------------------------------------
360 // We need at least 3 points to define a plane
361 if( (n = m_Points.Get_Count()) >= 3 && M.Create(n + 3, n + 3) && m_V.Create(n + 3) )
362 {
363 int i, j;
364 double a, b;
365 TSG_Point_Z Point;
366
367 //-------------------------------------------------
368 // Fill K (n x n, upper left of L) and calculate
369 // mean edge length from control points
370 //
371 // K is symmetrical so we really have to
372 // calculate only about half of the coefficients.
373 for(i=0, a=0.; i<n && (bSilent || SG_UI_Process_Set_Progress(i, n)); ++i )
374 {
375 Point = m_Points[i];
376
377 for(j=i+1; j<n; ++j)
378 {
379 b = _Get_hDistance(Point, m_Points[j]);
380 a += b * 2.; // same for upper & lower tri
381 M[i][j] = (M[j][i] = _Get_Base_Funtion(b));
382 }
383 }
384
385 a /= (double)(n*n);
386
387 //-------------------------------------------------
388 // Fill the rest of L
389 for(i=0; i<n; ++i)
390 {
391 // diagonal: reqularization parameters (lambda * a^2)
392 M[i][i] = Regularization * (a*a);
393
394 // P (n x 3, upper right)
395 M[i][n + 0] = 1.;
396 M[i][n + 1] = m_Points[i].x;
397 M[i][n + 2] = m_Points[i].y;
398
399 // P transposed (3 x n, bottom left)
400 M[n + 0][i] = 1.;
401 M[n + 1][i] = m_Points[i].x;
402 M[n + 2][i] = m_Points[i].y;
403 }
404
405 //-------------------------------------------------
406 // O (3 x 3, lower right)
407 for(i=n; i<n+3; ++i)
408 {
409 for(j=n; j<n+3; ++j)
410 {
411 M[i][j] = 0.;
412 }
413 }
414
415 //-------------------------------------------------
416 // Fill the right hand vector m_V
417 for(i=0; i<n; ++i)
418 {
419 m_V[i] = m_Points[i].z;
420 }
421
422 m_V[n + 0] = m_V[n + 1] = m_V[n + 2] = 0.;
423
424 //-------------------------------------------------
425 // Solve the linear system "inplace"
426 if( !bSilent )
427 {
428 SG_UI_Process_Set_Text(_TL("Thin Plate Spline: solving matrix"));
429 }
430
431 bResult = SG_Matrix_Solve(M, m_V, bSilent);
432 }
433
434 //-----------------------------------------------------
435 if( !bResult )
436 {
437 Destroy();
438 }
439
440 return( bResult );
441 }
442
443 //---------------------------------------------------------
Get_Value(double x,double y)444 double CSG_Thin_Plate_Spline::Get_Value(double x, double y)
445 {
446 if( m_V.Get_N() > 0 )
447 {
448 int n = m_Points.Get_Count();
449 double z = m_V[n + 0] + m_V[n + 1] * x + m_V[n + 2] * y;
450
451 for(int i=0; i<n; i++)
452 {
453 z += m_V[i] * _Get_Base_Funtion(m_Points[i], x, y);
454 }
455
456 return( z );
457 }
458
459 return( 0. );
460 }
461
462
463 ///////////////////////////////////////////////////////////
464 // //
465 // //
466 // //
467 ///////////////////////////////////////////////////////////
468
469 //---------------------------------------------------------
470