SAGA API  v9.5
mat_spline.cpp
Go to the documentation of this file.
1 
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 // //
51 
52 //---------------------------------------------------------
53 #include "mat_tools.h"
54 
55 
57 // //
58 // //
59 // //
61 
62 //---------------------------------------------------------
63 CSG_Spline::CSG_Spline(void)
64 {
65  m_bCreated = false;
66 }
67 
68 //---------------------------------------------------------
69 CSG_Spline::~CSG_Spline(void)
70 {
71  Destroy();
72 }
73 
74 
76 // //
78 
79 //---------------------------------------------------------
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 //---------------------------------------------------------
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 //---------------------------------------------------------
103 bool CSG_Spline::Create(double yA, double yB)
104 {
105  return( _Create(yA, yB) );
106 }
107 
108 
110 // //
112 
113 //---------------------------------------------------------
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 
124 // //
126 
127 //---------------------------------------------------------
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 
202 // //
204 
205 //---------------------------------------------------------
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 //---------------------------------------------------------
248 double CSG_Spline::Get_Value(double x)
249 {
250  Get_Value(x, x);
251 
252  return( x );
253 }
254 
255 
257 // //
258 // //
259 // //
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 
284 // //
285 // //
286 // //
288 
289 //---------------------------------------------------------
290 CSG_Thin_Plate_Spline::CSG_Thin_Plate_Spline(void)
291 {
292 }
293 
294 //---------------------------------------------------------
295 CSG_Thin_Plate_Spline::~CSG_Thin_Plate_Spline(void)
296 {
297  Destroy();
298 }
299 
300 
302 // //
304 
305 //---------------------------------------------------------
306 bool CSG_Thin_Plate_Spline::Destroy(void)
307 {
308  m_Points.Clear();
309  m_V.Destroy();
310 
311  return( true );
312 }
313 
314 
316 // //
318 
319 //---------------------------------------------------------
320 double CSG_Thin_Plate_Spline::_Get_hDistance(TSG_Point_3D A, TSG_Point_3D 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 //---------------------------------------------------------
329 double CSG_Thin_Plate_Spline::_Get_Base_Funtion(double x)
330 {
331  return( x > 0. ? x*x * log(x) : 0. );
332 }
333 
334 //---------------------------------------------------------
335 double CSG_Thin_Plate_Spline::_Get_Base_Funtion(TSG_Point_3D 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 
347 // //
349 
350 //---------------------------------------------------------
351 // Calculate Thin Plate Spline (TPS) weights from control points.
352 //
353 bool CSG_Thin_Plate_Spline::Create(double Regularization, bool bSilent)
354 {
355  sLong n = m_Points.Get_Count();
356 
357  if( n < 3 ) // We need at least 3 points to define a plane
358  {
359  Destroy();
360 
361  return( false );
362  }
363 
364  CSG_Matrix M;
365 
366  if( !M.Create(n + 3, n + 3) || !m_V.Create(n + 3) )
367  {
368  Destroy();
369 
370  return( false );
371  }
372 
373  //-----------------------------------------------------
374  // Fill K (n x n, upper left of L) and calculate
375  // mean edge length from control points
376  //
377  // K is symmetrical so we really have to
378  // calculate only about half of the coefficients.
379 
380  double a = 0.;
381 
382  for(sLong i=0; i<n && (bSilent || SG_UI_Process_Set_Progress(i, n)); ++i)
383  {
384  TSG_Point_3D Point = m_Points[i];
385 
386  for(sLong j=i+1; j<n; ++j)
387  {
388  double b = _Get_hDistance(Point, m_Points[j]); a += b * 2.; // same for upper & lower tri
389 
390  M[i][j] = (M[j][i] = _Get_Base_Funtion(b));
391  }
392  }
393 
394  a /= n*n;
395 
396  //-------------------------------------------------
397  // Fill the rest of L
398 
399  for(sLong i=0; i<n; ++i)
400  {
401  // P diagonal: reqularization parameters (lambda * a^2)
402  M[i][i] = Regularization * (a*a);
403 
404  // P (n x 3, upper right), P transposed (3 x n, bottom left)
405  M[i][n + 0] = M[n + 0][i] = 1.;
406  M[i][n + 1] = M[n + 1][i] = m_Points[i].x;
407  M[i][n + 2] = M[n + 2][i] = m_Points[i].y;
408 
409  // Fill the right hand vector m_V
410  m_V[i] = m_Points[i].z;
411  }
412 
413  for(sLong i=n; i<n+3; ++i)
414  {
415  for(sLong j=n; j<n+3; ++j)
416  {
417  M[i][j] = 0.; // O (3 x 3, lower right)
418  }
419  }
420 
421  m_V[n + 0] = m_V[n + 1] = m_V[n + 2] = 0.;
422 
423  //-------------------------------------------------
424  // Solve the linear system "inplace"
425 
426  if( !bSilent )
427  {
428  SG_UI_Process_Set_Text(_TL("Thin Plate Spline: solving matrix"));
429  }
430 
431  if( !SG_Matrix_Solve(M, m_V, bSilent) )
432  {
433  Destroy();
434 
435  return( false );
436  }
437 
438  return( true );
439 }
440 
441 //---------------------------------------------------------
442 double CSG_Thin_Plate_Spline::Get_Value(double x, double y)
443 {
444  if( m_V.Get_Size() > 0 )
445  {
446  sLong n = m_Points.Get_Count();
447  double z = m_V[n + 0] + m_V[n + 1] * x + m_V[n + 2] * y;
448 
449  for(sLong i=0; i<n; i++)
450  {
451  z += m_V[i] * _Get_Base_Funtion(m_Points[i], x, y);
452  }
453 
454  return( z );
455  }
456 
457  return( 0. );
458 }
459 
460 
462 // //
463 // //
464 // //
466 
467 //---------------------------------------------------------
CSG_Index
Definition: mat_tools.h:200
_TL
#define _TL(s)
Definition: api_core.h:1489
CSG_Vector::Get_Data
double * Get_Data(void) const
Definition: mat_tools.h:384
A
#define A
SSG_Point_3D
Definition: geo_tools.h:264
CSG_Points_3D::Get_Count
sLong Get_Count(void) const
Definition: geo_tools.h:336
CSG_Thin_Plate_Spline::Get_Value
double Get_Value(double x, double y)
Definition: mat_spline.cpp:442
CSG_Vector::Add_Row
bool Add_Row(double Value=0.)
Definition: mat_matrix.cpp:188
CSG_Points_3D::Clear
bool Clear(void)
Definition: geo_tools.h:326
CSG_Thin_Plate_Spline::Destroy
bool Destroy(void)
Definition: mat_spline.cpp:306
CSG_Spline::Create
bool Create(double *xValues, double *yValues, int nValues, double yA=1.0e30, double yB=1.0e30)
Definition: mat_spline.cpp:90
mat_tools.h
CSG_Spline::m_bCreated
bool m_bCreated
Definition: mat_tools.h:1438
CSG_Matrix::Create
bool Create(const CSG_Matrix &Matrix)
Definition: mat_matrix.cpp:836
CSG_Vector
Definition: mat_tools.h:360
CSG_Thin_Plate_Spline::Create
bool Create(double Regularization=0., bool bSilent=true)
Definition: mat_spline.cpp:353
CSG_Vector::Create
bool Create(const CSG_Vector &Vector)
Definition: mat_matrix.cpp:87
CSG_Spline::Destroy
void Destroy(void)
Definition: mat_spline.cpp:80
sLong
signed long long sLong
Definition: api_core.h:158
CSG_Spline::m_x
CSG_Vector m_x
Definition: mat_tools.h:1440
SG_UI_Process_Set_Text
void SG_UI_Process_Set_Text(const CSG_String &Text)
Definition: api_callback.cpp:324
CSG_Spline::m_y
CSG_Vector m_y
Definition: mat_tools.h:1440
CSG_Spline::_Create
bool _Create(double yA, double yB)
Definition: mat_spline.cpp:128
SG_Matrix_Solve
bool SG_Matrix_Solve(CSG_Matrix &Matrix, CSG_Vector &Vector, bool bSilent)
Definition: mat_matrix.cpp:1910
CSG_Spline::~CSG_Spline
virtual ~CSG_Spline(void)
Definition: mat_spline.cpp:69
CSG_Thin_Plate_Spline::CSG_Thin_Plate_Spline
CSG_Thin_Plate_Spline(void)
Definition: mat_spline.cpp:290
CSG_Vector::Get_Size
sLong Get_Size(void) const
Definition: mat_tools.h:380
CSG_Vector::Destroy
bool Destroy(void)
Definition: mat_matrix.cpp:130
CSG_Spline::Get_Value
bool Get_Value(double x, double &y)
Definition: mat_spline.cpp:206
B
#define B
SG_UI_Process_Set_Progress
bool SG_UI_Process_Set_Progress(int Position, int Range)
Definition: api_callback.cpp:256
CSG_Spline::CSG_Spline
CSG_Spline(void)
Definition: mat_spline.cpp:63
CSG_Spline::Add
void Add(double x, double y)
Definition: mat_spline.cpp:114
CSG_Spline::Get_Count
int Get_Count(void) const
Definition: mat_tools.h:1426
CSG_Matrix
Definition: mat_tools.h:478
CSG_Thin_Plate_Spline::~CSG_Thin_Plate_Spline
virtual ~CSG_Thin_Plate_Spline(void)
Definition: mat_spline.cpp:295
CSG_Spline::m_z
CSG_Vector m_z
Definition: mat_tools.h:1440