 |
SAGA API v2.0.8
|
|
SAGA API v2.0.8
|
00001 /********************************************************** 00002 * Version $Id: mat_spline.cpp 1173 2011-09-23 14:39:27Z oconrad $ 00003 *********************************************************/ 00004 00006 // // 00007 // SAGA // 00008 // // 00009 // System for Automated Geoscientific Analyses // 00010 // // 00011 // Application Programming Interface // 00012 // // 00013 // Library: SAGA_API // 00014 // // 00015 //-------------------------------------------------------// 00016 // // 00017 // mat_spline.cpp // 00018 // // 00019 // Copyright (C) 2005 by Olaf Conrad // 00020 // // 00021 //-------------------------------------------------------// 00022 // // 00023 // This file is part of 'SAGA - System for Automated // 00024 // Geoscientific Analyses'. // 00025 // // 00026 // This library is free software; you can redistribute // 00027 // it and/or modify it under the terms of the GNU Lesser // 00028 // General Public License as published by the Free // 00029 // Software Foundation, version 2.1 of the License. // 00030 // // 00031 // This library is distributed in the hope that it will // 00032 // be useful, but WITHOUT ANY WARRANTY; without even the // 00033 // implied warranty of MERCHANTABILITY or FITNESS FOR A // 00034 // PARTICULAR PURPOSE. See the GNU Lesser General Public // 00035 // License for more details. // 00036 // // 00037 // You should have received a copy of the GNU Lesser // 00038 // General Public License along with this program; if // 00039 // not, write to the Free Software Foundation, Inc., // 00040 // 59 Temple Place - Suite 330, Boston, MA 02111-1307, // 00041 // USA. // 00042 // // 00043 //-------------------------------------------------------// 00044 // // 00045 // contact: Olaf Conrad // 00046 // Institute of Geography // 00047 // University of Goettingen // 00048 // Goldschmidtstr. 5 // 00049 // 37077 Goettingen // 00050 // Germany // 00051 // // 00052 // e-mail: oconrad@saga-gis.org // 00053 // // 00055 00056 //--------------------------------------------------------- 00057 00058 00060 // // 00061 // // 00062 // // 00064 00065 //--------------------------------------------------------- 00066 #include "mat_tools.h" 00067 00068 00070 // // 00071 // // 00072 // // 00074 00075 //--------------------------------------------------------- 00076 CSG_Spline::CSG_Spline(void) 00077 { 00078 m_Values = NULL; 00079 m_nValues = 0; 00080 m_nBuffer = 0; 00081 m_bCreated = false; 00082 } 00083 00084 //--------------------------------------------------------- 00085 CSG_Spline::~CSG_Spline(void) 00086 { 00087 Destroy(); 00088 } 00089 00090 00092 // // 00094 00095 //--------------------------------------------------------- 00096 void CSG_Spline::Destroy(void) 00097 { 00098 if( m_Values ) 00099 { 00100 SG_Free(m_Values); 00101 00102 m_Values = NULL; 00103 m_nValues = 0; 00104 m_nBuffer = 0; 00105 m_bCreated = false; 00106 } 00107 } 00108 00109 //--------------------------------------------------------- 00110 bool CSG_Spline::Create(double *xValues, double *yValues, int nValues, double yA, double yB) 00111 { 00112 Destroy(); 00113 00114 for(int i=0; i<nValues; i++) 00115 { 00116 Add(xValues[i], yValues[i]); 00117 } 00118 00119 return( _Create(yA, yB) ); 00120 } 00121 00122 //--------------------------------------------------------- 00123 bool CSG_Spline::Create(double yA, double yB) 00124 { 00125 return( _Create(yA, yB) ); 00126 } 00127 00128 00130 // // 00132 00133 //--------------------------------------------------------- 00134 void CSG_Spline::Add(double x, double y) 00135 { 00136 m_bCreated = false; 00137 00138 //----------------------------------------------------- 00139 if( m_nValues >= m_nBuffer ) 00140 { 00141 m_nBuffer += 64; 00142 m_Values = (TSG_Point_Z *)SG_Realloc(m_Values, m_nBuffer * sizeof(TSG_Point_Z)); 00143 } 00144 00145 m_nValues++; 00146 00147 //----------------------------------------------------- 00148 if( m_nValues == 1 ) 00149 { 00150 m_Values[0].x = x; 00151 m_Values[0].y = y; 00152 } 00153 else 00154 { 00155 int i, iAdd; 00156 00157 for(iAdd=0; iAdd<m_nValues-1 && m_Values[iAdd].x<x; iAdd++) {} 00158 00159 for(i=m_nValues-1; i>iAdd; i--) 00160 { 00161 m_Values[i] = m_Values[i - 1]; 00162 } 00163 00164 m_Values[iAdd].x = x; 00165 m_Values[iAdd].y = y; 00166 } 00167 } 00168 00169 00171 // // 00173 00174 //--------------------------------------------------------- 00175 bool CSG_Spline::_Create(double yA, double yB) 00176 { 00177 int i, k; 00178 double p, qn, sig, un, *u; 00179 00180 if( m_nValues > 2 ) 00181 { 00182 m_bCreated = true; 00183 u = (double *)SG_Malloc(m_nValues * sizeof(double)); 00184 00185 if( yA > 0.99e30 ) 00186 { 00187 m_Values[0].z = u[0] = 0.0; 00188 } 00189 else 00190 { 00191 m_Values[0].z = -0.5; 00192 u[0] = (3.0 / (m_Values[1].x - m_Values[0].x)) 00193 * ((m_Values[1].y - m_Values[0].y) / (m_Values[1].x - m_Values[0].x) - yA); 00194 } 00195 00196 for(i=1; i<m_nValues-1; i++) 00197 { 00198 sig = (m_Values[i].x - m_Values[i - 1].x) / (m_Values[i + 1].x - m_Values[i - 1].x); 00199 p = sig * m_Values[i - 1].z + 2.0; 00200 m_Values[i].z = (sig - 1.0) / p; 00201 u[i] = (m_Values[i + 1].y - m_Values[i ].y) / (m_Values[i + 1].x - m_Values[i ].x) 00202 - (m_Values[i ].y - m_Values[i - 1].y) / (m_Values[i ].x - m_Values[i - 1].x); 00203 u[i] = (6.0 * u[i] / (m_Values[i + 1].x - m_Values[i - 1].x) - sig * u[i - 1]) / p; 00204 } 00205 00206 if( yB > 0.99e30 ) 00207 { 00208 qn = un = 0.0; 00209 } 00210 else 00211 { 00212 qn = 0.5; 00213 un = (3.0 / (m_Values[m_nValues - 1].x - m_Values[m_nValues - 2].x)) 00214 * (yB - (m_Values[m_nValues - 1].y - m_Values[m_nValues - 2].y) 00215 / (m_Values[m_nValues - 1].x - m_Values[m_nValues - 2].x)); 00216 } 00217 00218 m_Values[m_nValues - 1].z = (un - qn * u[m_nValues - 2]) / (qn * m_Values[m_nValues - 2].z + 1.0); 00219 00220 for(k=m_nValues-2; k>=0; k--) 00221 { 00222 m_Values[k].z = m_Values[k].z * m_Values[k + 1].z + u[k]; 00223 } 00224 00225 SG_Free(u); 00226 00227 return( true ); 00228 } 00229 00230 return( false ); 00231 } 00232 00233 00235 // // 00237 00238 //--------------------------------------------------------- 00239 bool CSG_Spline::Get_Value(double x, double &y) 00240 { 00241 if( m_bCreated || Create() ) 00242 { 00243 int klo, khi, k; 00244 double h, b, a; 00245 00246 klo = 0; 00247 khi = m_nValues - 1; 00248 00249 while( khi - klo > 1 ) 00250 { 00251 k = (khi+klo) >> 1; 00252 00253 if( m_Values[k].x > x ) 00254 { 00255 khi = k; 00256 } 00257 else 00258 { 00259 klo = k; 00260 } 00261 } 00262 00263 h = m_Values[khi].x - m_Values[klo].x; 00264 00265 if( h != 0.0 ) 00266 { 00267 a = (m_Values[khi].x - x) / h; 00268 b = (x - m_Values[klo].x) / h; 00269 00270 y = a * m_Values[klo].y + b * m_Values[khi].y 00271 + ((a*a*a - a) * m_Values[klo].z + (b*b*b - b) * m_Values[khi].z) * (h*h) / 6.0; 00272 00273 return( true ); 00274 } 00275 } 00276 00277 return( false ); 00278 } 00279 00280 //--------------------------------------------------------- 00281 double CSG_Spline::Get_Value(double x) 00282 { 00283 Get_Value(x, x); 00284 00285 return( x ); 00286 } 00287 00288 00290 // // 00291 // // 00292 // // 00294 00295 //--------------------------------------------------------- 00296 // 00297 // Based on: 00298 // 00299 // Thin Plate Spline demo/example in C++ 00300 // Copyright (C) 2003, 2005 by Jarno Elonen 00301 // 00302 // Permission to use, copy, modify, distribute and sell this software 00303 // and its documentation for any purpose is hereby granted without fee, 00304 // provided that the above copyright notice appear in all copies and 00305 // that both that copyright notice and this permission notice appear 00306 // in supporting documentation. The authors make no representations 00307 // about the suitability of this software for any purpose. 00308 // It is provided "as is" without express or implied warranty. 00309 // 00310 // Reference: 00311 // - Donato, G., Belongie, S. (2002): 00312 // 'Approximation Methods for Thin Plate Spline Mappings and Principal Warps' 00313 // 00314 //--------------------------------------------------------- 00315 00317 // // 00318 // // 00319 // // 00321 00322 //--------------------------------------------------------- 00323 CSG_Thin_Plate_Spline::CSG_Thin_Plate_Spline(void) 00324 { 00325 } 00326 00327 //--------------------------------------------------------- 00328 CSG_Thin_Plate_Spline::~CSG_Thin_Plate_Spline(void) 00329 { 00330 Destroy(); 00331 } 00332 00333 00335 // // 00337 00338 //--------------------------------------------------------- 00339 bool CSG_Thin_Plate_Spline::Destroy(void) 00340 { 00341 m_Points.Clear(); 00342 m_V.Destroy(); 00343 00344 return( true ); 00345 } 00346 00347 00349 // // 00351 00352 //--------------------------------------------------------- 00353 double CSG_Thin_Plate_Spline::_Get_hDistance(TSG_Point_Z A, TSG_Point_Z B) 00354 { 00355 A.x -= B.x; 00356 A.y -= B.y; 00357 00358 return( sqrt(A.x*A.x + A.y*A.y) ); 00359 } 00360 00361 //--------------------------------------------------------- 00362 double CSG_Thin_Plate_Spline::_Get_Base_Funtion(double x) 00363 { 00364 return( x > 0.0 ? x*x * log(x) : 0.0 ); 00365 } 00366 00367 //--------------------------------------------------------- 00368 double CSG_Thin_Plate_Spline::_Get_Base_Funtion(TSG_Point_Z A, double x, double y) 00369 { 00370 x -= A.x; 00371 y -= A.y; 00372 x = sqrt(x*x + y*y); 00373 00374 return( x > 0.0 ? x*x * log(x) : 0.0 ); 00375 } 00376 00377 00379 // // 00381 00382 //--------------------------------------------------------- 00383 // Calculate Thin Plate Spline (TPS) weights from control points. 00384 // 00385 bool CSG_Thin_Plate_Spline::Create(double Regularization, bool bSilent) 00386 { 00387 bool bResult = false; 00388 int n; 00389 00390 //----------------------------------------------------- 00391 // You We need at least 3 points to define a plane 00392 if( (n = m_Points.Get_Count()) >= 3 ) 00393 { 00394 int i, j; 00395 double a, b; 00396 TSG_Point_Z Point; 00397 CSG_Matrix M; 00398 00399 //------------------------------------------------- 00400 // Allocate the matrix and vector 00401 M .Create(n + 3, n + 3); 00402 m_V .Create(n + 3); 00403 00404 //------------------------------------------------- 00405 // Fill K (n x n, upper left of L) and calculate 00406 // mean edge length from control points 00407 // 00408 // K is symmetrical so we really have to 00409 // calculate only about half of the coefficients. 00410 for(i=0, a=0.0; i<n && (bSilent || SG_UI_Process_Set_Progress(i, n)); ++i ) 00411 { 00412 Point = m_Points[i]; 00413 00414 for(j=i+1; j<n; ++j) 00415 { 00416 b = _Get_hDistance(Point, m_Points[j]); 00417 a += b * 2.0; // same for upper & lower tri 00418 M[i][j] = (M[j][i] = _Get_Base_Funtion(b)); 00419 } 00420 } 00421 00422 a /= (double)(n*n); 00423 00424 //------------------------------------------------- 00425 // Fill the rest of L 00426 for(i=0; i<n; ++i) 00427 { 00428 // diagonal: reqularization parameters (lambda * a^2) 00429 M[i][i] = Regularization * (a*a); 00430 00431 // P (n x 3, upper right) 00432 M[i][n + 0] = 1.0; 00433 M[i][n + 1] = m_Points[i].x; 00434 M[i][n + 2] = m_Points[i].y; 00435 00436 // P transposed (3 x n, bottom left) 00437 M[n + 0][i] = 1.0; 00438 M[n + 1][i] = m_Points[i].x; 00439 M[n + 2][i] = m_Points[i].y; 00440 } 00441 00442 //------------------------------------------------- 00443 // O (3 x 3, lower right) 00444 for(i=n; i<n+3; ++i) 00445 { 00446 for(j=n; j<n+3; ++j) 00447 { 00448 M[i][j] = 0.0; 00449 } 00450 } 00451 00452 //------------------------------------------------- 00453 // Fill the right hand vector m_V 00454 for(i=0; i<n; ++i) 00455 { 00456 m_V[i] = m_Points[i].z; 00457 } 00458 00459 m_V[n + 0] = m_V[n + 1] = m_V[n + 2] = 0.0; 00460 00461 //------------------------------------------------- 00462 // Solve the linear system "inplace" 00463 if( !bSilent ) 00464 { 00465 SG_UI_Process_Set_Text(LNG("Thin Plate Spline: solving matrix")); 00466 } 00467 00468 bResult = SG_Matrix_Solve(M, m_V, bSilent); 00469 } 00470 00471 //----------------------------------------------------- 00472 if( !bResult ) 00473 { 00474 Destroy(); 00475 } 00476 00477 return( bResult ); 00478 } 00479 00480 //--------------------------------------------------------- 00481 double CSG_Thin_Plate_Spline::Get_Value(double x, double y) 00482 { 00483 if( m_V.Get_N() > 0 ) 00484 { 00485 int n = m_Points.Get_Count(); 00486 double z = m_V[n + 0] + m_V[n + 1] * x + m_V[n + 2] * y; 00487 00488 for(int i=0; i<n; i++) 00489 { 00490 z += m_V[i] * _Get_Base_Funtion(m_Points[i], x, y); 00491 } 00492 00493 return( z ); 00494 } 00495 00496 return( 0.0 ); 00497 } 00498 00499 00501 // // 00502 // // 00503 // // 00505 00506 //---------------------------------------------------------
1.7.4