1
2 ///////////////////////////////////////////////////////////
3 // //
4 // SAGA //
5 // //
6 // System for Automated Geoscientific Analyses //
7 // //
8 // Tool Library //
9 // grid_spline //
10 // //
11 //-------------------------------------------------------//
12 // //
13 // Gridding_Spline_MBA.cpp //
14 // //
15 // Copyright (C) 2006 by //
16 // Olaf Conrad //
17 // //
18 //-------------------------------------------------------//
19 // //
20 // This file is part of 'SAGA - System for Automated //
21 // Geoscientific Analyses'. SAGA is free software; you //
22 // can redistribute it and/or modify it under the terms //
23 // of the GNU General Public License as published by the //
24 // Free Software Foundation, either version 2 of the //
25 // License, or (at your option) any later version. //
26 // //
27 // SAGA is distributed in the hope that it will be //
28 // useful, but WITHOUT ANY WARRANTY; without even the //
29 // implied warranty of MERCHANTABILITY or FITNESS FOR A //
30 // PARTICULAR PURPOSE. See the GNU General Public //
31 // License for more details. //
32 // //
33 // You should have received a copy of the GNU General //
34 // Public License along with this program; if not, see //
35 // <http://www.gnu.org/licenses/>. //
36 // //
37 //-------------------------------------------------------//
38 // //
39 // e-mail: oconrad@saga-gis.org //
40 // //
41 // contact: Olaf Conrad //
42 // Institute of Geography //
43 // University of Goettingen //
44 // Goldschmidtstr. 5 //
45 // 37077 Goettingen //
46 // Germany //
47 // //
48 ///////////////////////////////////////////////////////////
49
50 //---------------------------------------------------------
51 #include "Gridding_Spline_MBA.h"
52
53
54 ///////////////////////////////////////////////////////////
55 // //
56 // //
57 // //
58 ///////////////////////////////////////////////////////////
59
60 //---------------------------------------------------------
CGridding_Spline_MBA(void)61 CGridding_Spline_MBA::CGridding_Spline_MBA(void)
62 : CGridding_Spline_Base()
63 {
64 Set_Name (_TL("Multilevel B-Spline"));
65
66 Set_Author ("O.Conrad (c) 2006");
67
68 Set_Description (_TW(
69 "Multilevel B-spline algorithm for spatial interpolation of scattered data "
70 "as proposed by Lee, Wolberg and Shin (1997).\n"
71 "The algorithm makes use of a coarse-to-fine hierarchy of control lattices to "
72 "generate a sequence of bicubic B-spline functions, whose sum approaches the "
73 "desired interpolation function. Performance gains are realized by using "
74 "B-spline refinement to reduce the sum of these functions into one equivalent "
75 "B-spline function. "
76 "\n\n"
77 "The 'Maximum Level' determines the maximum size of the final B-spline matrix "
78 "and increases exponential with each level. Where level=10 requires about 1mb "
79 "level=12 needs about 16mb and level=14 about 256mb(!) of additional memory. "
80 ));
81
82 Add_Reference(
83 "Lee, S., Wolberg, G., Shin, S.Y.", "1997",
84 "Scattered Data Interpolation with Multilevel B-Splines",
85 "IEEE Transactions On Visualisation And Computer Graphics, Vol.3, No.3., p.228-244.",
86 SG_T("https://www.researchgate.net/profile/George_Wolberg/publication/3410822_Scattered_Data_Interpolation_with_Multilevel_B-Splines/links/00b49518719ac9f08a000000/Scattered-Data-Interpolation-with-Multilevel-B-Splines.pdf"),
87 SG_T("ResearchGate")
88 );
89
90 //-----------------------------------------------------
91 Parameters.Add_Choice(
92 "", "METHOD" , _TL("Refinement"),
93 _TL(""),
94 CSG_String::Format("%s|%s",
95 _TL("no"),
96 _TL("yes")
97 ), 0
98 );
99
100 Parameters.Add_Double(
101 "", "EPSILON" , _TL("Threshold Error"),
102 _TL(""),
103 0.0001, 0.0, true
104 );
105
106 Parameters.Add_Int(
107 "", "LEVEL_MAX" , _TL("Maximum Level"),
108 _TL(""),
109 11, 1, true, 14, true
110 );
111
112 Parameters.Add_Bool(
113 "", "UPDATE" , _TL("Update View"),
114 _TL(""),
115 false
116 )->Set_UseInCMD(false);
117 }
118
119
120 ///////////////////////////////////////////////////////////
121 // //
122 ///////////////////////////////////////////////////////////
123
124 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)125 int CGridding_Spline_MBA::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
126 {
127 if( pParameter->Cmp_Identifier("METHOD") )
128 {
129 pParameters->Set_Enabled("UPDATE", pParameter->asInt() == 0); // no performance gain with refinement!
130 }
131
132 return( CGridding_Spline_Base::On_Parameters_Enable(pParameters, pParameter) );
133 }
134
135
136 ///////////////////////////////////////////////////////////
137 // //
138 ///////////////////////////////////////////////////////////
139
140 //---------------------------------------------------------
On_Execute(void)141 bool CGridding_Spline_MBA::On_Execute(void)
142 {
143 bool bResult = false;
144
145 if( Initialize(m_Points, true, true) )
146 {
147 m_Epsilon = Parameters("EPSILON")->asDouble();
148
149 double Cellsize = M_GET_MAX(m_pGrid->Get_XRange(), m_pGrid->Get_YRange());
150
151 switch( Parameters("METHOD")->asInt() )
152 {
153 case 0: bResult = _Set_MBA (Cellsize); break;
154 default: bResult = _Set_MBA_Refinement(Cellsize); break;
155 }
156
157 m_Points.Clear();
158
159 Finalize(true);
160 }
161
162 return( bResult );
163 }
164
165
166 ///////////////////////////////////////////////////////////
167 // //
168 ///////////////////////////////////////////////////////////
169
170 //---------------------------------------------------------
_Set_MBA(double Cellsize)171 bool CGridding_Spline_MBA::_Set_MBA(double Cellsize)
172 {
173 CSG_Grid Phi;
174
175 bool bContinue = true;
176
177 int Levels = Parameters("LEVEL_MAX")->asInt();
178
179 for(int Level=0; bContinue && Level<Levels && Process_Get_Okay(false); Level++, Cellsize/=2.)
180 {
181 bContinue = BA_Set_Phi(Phi, Cellsize) && _Get_Difference(Phi, Level);
182
183 BA_Set_Grid(Phi, Level > 0);
184
185 if( Parameters("UPDATE")->asBool() )
186 {
187 DataObject_Update(m_pGrid, true);
188 }
189 }
190
191 return( true );
192 }
193
194
195 ///////////////////////////////////////////////////////////
196 // //
197 ///////////////////////////////////////////////////////////
198
199 //---------------------------------------------------------
_Set_MBA_Refinement(double Cellsize)200 bool CGridding_Spline_MBA::_Set_MBA_Refinement(double Cellsize)
201 {
202 CSG_Grid Phi[2];
203
204 bool bContinue = true;
205
206 int Levels = Parameters("LEVEL_MAX")->asInt(), i = 0;
207
208 for(int Level=0; bContinue && Level<Levels && Process_Get_Okay(false); Level++, Cellsize/=2.)
209 {
210 i = Level % 2;
211
212 bContinue = BA_Set_Phi(Phi[i], Cellsize) && _Get_Difference(Phi[i], Level);
213
214 _Set_MBA_Refinement(Phi[(i + 1) % 2], Phi[i]);
215 }
216
217 BA_Set_Grid(Phi[i]);
218
219 return( true );
220 }
221
222 //---------------------------------------------------------
_Set_MBA_Refinement(const CSG_Grid & Psi_0,CSG_Grid & Psi_1)223 bool CGridding_Spline_MBA::_Set_MBA_Refinement(const CSG_Grid &Psi_0, CSG_Grid &Psi_1)
224 {
225 if( 2 * (Psi_0.Get_NX() - 4) != (Psi_1.Get_NX() - 4)
226 || 2 * (Psi_0.Get_NY() - 4) != (Psi_1.Get_NY() - 4) )
227 {
228 return( false );
229 }
230
231 #pragma omp parallel for
232 for(int y=0; y<Psi_0.Get_NY(); y++)
233 {
234 int yy = 2 * y - 1;
235
236 for(int x=0, xx=-1; x<Psi_0.Get_NX(); x++, xx+=2)
237 {
238 double a[3][3];
239
240 for(int iy=0, jy=y-1; iy<3; iy++, jy++)
241 {
242 for(int ix=0, jx=x-1; ix<3; ix++, jx++)
243 {
244 a[ix][iy] = Psi_0.is_InGrid(jx, jy, false) ? Psi_0.asDouble(jx, jy) : 0.0;
245 }
246 }
247
248 #define SET_PSI(x, y, z) if( Psi_1.is_InGrid(x, y) ) { Psi_1.Add_Value(x, y, z); }
249
250 SET_PSI(xx + 0, yy + 0,
251 ( a[0][0] + a[0][2] + a[2][0] + a[2][2]
252 + (a[0][1] + a[1][0] + a[1][2] + a[2][1]) * 6.0
253 + a[1][1] * 36.0
254 ) / 64.0
255 );
256
257 SET_PSI(xx + 0, yy + 1,
258 ( a[0][1] + a[0][2] + a[2][1] + a[2][2]
259 + (a[1][1] + a[1][2]) * 6.0
260 ) / 16.0
261 );
262
263 SET_PSI(xx + 1, yy + 0,
264 ( a[1][0] + a[1][2] + a[2][0] + a[2][2]
265 + (a[1][1] + a[2][1]) * 6.0
266 ) / 16.0
267 );
268
269 SET_PSI(xx + 1, yy + 1,
270 ( a[1][1] + a[1][2] + a[2][1] + a[2][2]
271 ) / 4.0
272 );
273 }
274 }
275
276 return( true );
277 }
278
279
280 ///////////////////////////////////////////////////////////
281 // //
282 ///////////////////////////////////////////////////////////
283
284 //---------------------------------------------------------
_Get_Difference(const CSG_Grid & Phi,int Level)285 bool CGridding_Spline_MBA::_Get_Difference(const CSG_Grid &Phi, int Level)
286 {
287 CSG_Simple_Statistics Differences;
288
289 for(int i=0; i<m_Points.Get_Count(); i++)
290 {
291 TSG_Point_Z p = m_Points[i];
292
293 p.x = (p.x - Phi.Get_XMin()) / Phi.Get_Cellsize();
294 p.y = (p.y - Phi.Get_YMin()) / Phi.Get_Cellsize();
295 p.z = p.z - BA_Get_Phi(Phi, p.x, p.y);
296
297 m_Points[i].z = p.z;
298
299 if( fabs(p.z) > m_Epsilon )
300 {
301 Differences += fabs(p.z);
302 }
303 }
304
305 //-----------------------------------------------------
306 Message_Fmt("\n%s:%d %s:%d %s:%f %s:%f",
307 _TL("level" ), Level + 1,
308 _TL("errors" ), (int)Differences.Get_Count (),
309 _TL("maximum"), Differences.Get_Maximum(),
310 _TL("mean" ), Differences.Get_Mean ()
311 );
312
313 Process_Set_Text(CSG_String::Format("%s %d [%d]", _TL("Level"), Level + 1, (int)Differences.Get_Count()));
314
315 return( Differences.Get_Maximum() > m_Epsilon );
316 }
317
318
319 ///////////////////////////////////////////////////////////
320 // //
321 ///////////////////////////////////////////////////////////
322
323 //---------------------------------------------------------
BA_Get_B(int i,double d) const324 inline double CGridding_Spline_MBA::BA_Get_B(int i, double d) const
325 {
326 switch( i )
327 {
328 case 0: d = 1. - d; return( d*d*d / 6. );
329
330 case 1: return( ( 3. * d*d*d - 6. * d*d + 4.) / 6. );
331
332 case 2: return( (-3. * d*d*d + 3. * d*d + 3. * d + 1.) / 6. );
333
334 case 3: return( d*d*d / 6. );
335
336 default: return( 0. );
337 }
338 }
339
340 //---------------------------------------------------------
BA_Set_Phi(CSG_Grid & Phi,double Cellsize)341 bool CGridding_Spline_MBA::BA_Set_Phi(CSG_Grid &Phi, double Cellsize)
342 {
343 int n = 4 + (int)(M_GET_MAX(m_pGrid->Get_XRange(), m_pGrid->Get_YRange()) / Cellsize);
344
345 Phi.Create(SG_DATATYPE_Float, n, n, Cellsize, m_pGrid->Get_XMin(), m_pGrid->Get_YMin());
346
347 CSG_Grid Delta(Phi.Get_System());
348
349 //-----------------------------------------------------
350 for(int i=0; i<m_Points.Get_Count(); i++)
351 {
352 TSG_Point_Z p = m_Points[i];
353
354 int x = (int)(p.x = (p.x - Phi.Get_XMin()) / Phi.Get_Cellsize());
355 int y = (int)(p.y = (p.y - Phi.Get_YMin()) / Phi.Get_Cellsize());
356
357 if( x >= 0 && x < Phi.Get_NX() - 3 && y >= 0 && y < Phi.Get_NY() - 3 )
358 {
359 int iy; double W[4][4], SW2 = 0.0;
360
361 for(iy=0; iy<4; iy++) // compute W[k,l] and Sum[a=0-3, b=0-3](W�[a,b])
362 {
363 double wy = BA_Get_B(iy, p.y - y);
364
365 for(int ix=0; ix<4; ix++)
366 {
367 SW2 += SG_Get_Square(W[iy][ix] = wy * BA_Get_B(ix, p.x - x));
368 }
369 }
370
371 if( SW2 > 0.0 )
372 {
373 p.z /= SW2;
374
375 for(iy=0; iy<4; iy++)
376 {
377 for(int ix=0; ix<4; ix++)
378 {
379 double wxy = W[iy][ix];
380
381 Delta.Add_Value(x + ix, y + iy, wxy*wxy*wxy * p.z); // numerator
382 Phi .Add_Value(x + ix, y + iy, wxy*wxy ); // denominator
383 }
384 }
385 }
386 }
387 }
388
389 //-----------------------------------------------------
390 #pragma omp parallel for
391 for(int y=0; y<Phi.Get_NY(); y++)
392 {
393 for(int x=0; x<Phi.Get_NX(); x++)
394 {
395 double z = Phi.asDouble(x, y);
396
397 if( z != 0. )
398 {
399 Phi.Set_Value(x, y, Delta.asDouble(x, y) / z);
400 }
401 }
402 }
403
404 //-----------------------------------------------------
405 return( true );
406 }
407
408 //---------------------------------------------------------
BA_Get_Phi(const CSG_Grid & Phi,double px,double py) const409 double CGridding_Spline_MBA::BA_Get_Phi(const CSG_Grid &Phi, double px, double py) const
410 {
411 double z = 0.0;
412
413 int x = (int)px; px -= x;
414 int y = (int)py; py -= y;
415
416 if( x >= 0 && x < Phi.Get_NX() - 3 && y >= 0 && y < Phi.Get_NY() - 3 )
417 {
418 for(int iy=0; iy<4; iy++)
419 {
420 double by = BA_Get_B(iy, py);
421
422 for(int ix=0; ix<4; ix++)
423 {
424 z += by * BA_Get_B(ix, px) * Phi.asDouble(x + ix, y + iy);
425 }
426 }
427 }
428
429 return( z );
430 }
431
432 //---------------------------------------------------------
BA_Set_Grid(const CSG_Grid & Phi,bool bAdd)433 void CGridding_Spline_MBA::BA_Set_Grid(const CSG_Grid &Phi, bool bAdd)
434 {
435 double d = m_pGrid->Get_Cellsize() / Phi.Get_Cellsize();
436
437 #pragma omp parallel for
438 for(int y=0; y<m_pGrid->Get_NY(); y++)
439 {
440 double py = d * y;
441
442 for(int x=0; x<m_pGrid->Get_NX(); x++)
443 {
444 double px = d * x;
445
446 if( bAdd )
447 { m_pGrid->Add_Value(x, y, BA_Get_Phi(Phi, px, py)); }
448 else
449 { m_pGrid->Set_Value(x, y, BA_Get_Phi(Phi, px, py)); }
450 }
451 }
452 }
453
454
455 ///////////////////////////////////////////////////////////
456 // //
457 // //
458 // //
459 ///////////////////////////////////////////////////////////
460
461 //---------------------------------------------------------
462