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_Grid.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_Grid.h"
52
53
54 ///////////////////////////////////////////////////////////
55 // //
56 // //
57 // //
58 ///////////////////////////////////////////////////////////
59
60 //---------------------------------------------------------
CGridding_Spline_MBA_Grid(void)61 CGridding_Spline_MBA_Grid::CGridding_Spline_MBA_Grid(void)
62 : CGridding_Spline_Base(true)
63 {
64 Set_Name (_TL("Multilevel B-Spline from Grid Points"));
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). "
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. Large 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., 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 Parameters.Add_Choice("TARGET",
119 "DATATYPE" , _TL("Data Type"),
120 _TL(""),
121 CSG_String::Format("%s|%s",
122 _TL("same as input grid"),
123 SG_Data_Type_Get_Name(SG_DATATYPE_Float ).c_str()
124 ), 0
125 );
126 }
127
128
129 ///////////////////////////////////////////////////////////
130 // //
131 ///////////////////////////////////////////////////////////
132
133 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)134 int CGridding_Spline_MBA_Grid::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
135 {
136 if( pParameter->Cmp_Identifier("METHOD") )
137 {
138 pParameters->Set_Enabled("UPDATE", pParameter->asInt() == 0); // no performance gain with refinement!
139 }
140
141 return( CGridding_Spline_Base::On_Parameters_Enable(pParameters, pParameter) );
142 }
143
144
145 ///////////////////////////////////////////////////////////
146 // //
147 ///////////////////////////////////////////////////////////
148
149 //---------------------------------------------------------
On_Execute(void)150 bool CGridding_Spline_MBA_Grid::On_Execute(void)
151 {
152 bool bResult = false;
153
154 if( Initialize() )
155 {
156 if( Parameters("DATATYPE")->asInt() == 0 )
157 {
158 m_Points.Create(*Parameters("GRID")->asGrid());
159 }
160 else
161 {
162 m_Points.Create(Parameters("GRID")->asGrid(), SG_DATATYPE_Float);
163 m_Points.Assign(Parameters("GRID")->asGrid());
164 }
165
166 m_Points.Add(-Parameters("GRID")->asGrid()->Get_Mean()); // detrend
167
168 m_Epsilon = Parameters("EPSILON")->asDouble();
169
170 double Cellsize = M_GET_MAX(m_pGrid->Get_XRange(), m_pGrid->Get_YRange());
171
172 switch( Parameters("METHOD")->asInt() )
173 {
174 case 0: bResult = _Set_MBA (Cellsize); break;
175 default: bResult = _Set_MBA_Refinement(Cellsize); break;
176 }
177
178 m_Points.Destroy();
179
180 Finalize(true);
181 }
182
183 return( bResult );
184 }
185
186
187 ///////////////////////////////////////////////////////////
188 // //
189 ///////////////////////////////////////////////////////////
190
191 //---------------------------------------------------------
_Set_MBA(double Cellsize)192 bool CGridding_Spline_MBA_Grid::_Set_MBA(double Cellsize)
193 {
194 CSG_Grid Phi;
195
196 bool bContinue = true;
197
198 int Levels = Parameters("LEVEL_MAX")->asInt();
199
200 for(int Level=0; bContinue && Level<Levels && Process_Get_Okay(false); Level++, Cellsize/=2.)
201 {
202 bContinue = BA_Set_Phi(Phi, Cellsize) && _Get_Difference(Phi, Level);
203
204 BA_Set_Grid(Phi, Level > 0);
205
206 if( Parameters("UPDATE")->asBool() )
207 {
208 DataObject_Update(m_pGrid, true);
209 }
210 }
211
212 return( true );
213 }
214
215
216 ///////////////////////////////////////////////////////////
217 // //
218 ///////////////////////////////////////////////////////////
219
220 //---------------------------------------------------------
_Set_MBA_Refinement(double Cellsize)221 bool CGridding_Spline_MBA_Grid::_Set_MBA_Refinement(double Cellsize)
222 {
223 CSG_Grid Phi[2];
224
225 bool bContinue = true;
226
227 int Levels = Parameters("LEVEL_MAX")->asInt(), i = 0;
228
229 for(int Level=0; bContinue && Level<Levels && Process_Get_Okay(false); Level++, Cellsize/=2.)
230 {
231 i = Level % 2;
232
233 bContinue = BA_Set_Phi(Phi[i], Cellsize) && _Get_Difference(Phi[i], Level);
234
235 _Set_MBA_Refinement(Phi[(i + 1) % 2], Phi[i]);
236 }
237
238 BA_Set_Grid(Phi[i]);
239
240 return( true );
241 }
242
243 //---------------------------------------------------------
_Set_MBA_Refinement(const CSG_Grid & Psi_0,CSG_Grid & Psi_1)244 bool CGridding_Spline_MBA_Grid::_Set_MBA_Refinement(const CSG_Grid &Psi_0, CSG_Grid &Psi_1)
245 {
246 if( 2 * (Psi_0.Get_NX() - 4) != (Psi_1.Get_NX() - 4)
247 || 2 * (Psi_0.Get_NY() - 4) != (Psi_1.Get_NY() - 4) )
248 {
249 return( false );
250 }
251
252 #pragma omp parallel for
253 for(int y=0; y<Psi_0.Get_NY(); y++)
254 {
255 int yy = 2 * y - 1;
256
257 for(int x=0, xx=-1; x<Psi_0.Get_NX(); x++, xx+=2)
258 {
259 double a[3][3];
260
261 for(int iy=0, jy=y-1; iy<3; iy++, jy++)
262 {
263 for(int ix=0, jx=x-1; ix<3; ix++, jx++)
264 {
265 a[ix][iy] = Psi_0.is_InGrid(jx, jy, false) ? Psi_0.asDouble(jx, jy) : 0.;
266 }
267 }
268
269 #define SET_PSI(x, y, z) if( Psi_1.is_InGrid(x, y) ) { Psi_1.Add_Value(x, y, z); }
270
271 SET_PSI(xx + 0, yy + 0,
272 ( a[0][0] + a[0][2] + a[2][0] + a[2][2]
273 + (a[0][1] + a[1][0] + a[1][2] + a[2][1]) * 6.
274 + a[1][1] * 36.
275 ) / 64.
276 );
277
278 SET_PSI(xx + 0, yy + 1,
279 ( a[0][1] + a[0][2] + a[2][1] + a[2][2]
280 + (a[1][1] + a[1][2]) * 6.
281 ) / 16.
282 );
283
284 SET_PSI(xx + 1, yy + 0,
285 ( a[1][0] + a[1][2] + a[2][0] + a[2][2]
286 + (a[1][1] + a[2][1]) * 6.
287 ) / 16.
288 );
289
290 SET_PSI(xx + 1, yy + 1,
291 ( a[1][1] + a[1][2] + a[2][1] + a[2][2]
292 ) / 4.
293 );
294 }
295 }
296
297 return( true );
298 }
299
300
301 ///////////////////////////////////////////////////////////
302 // //
303 ///////////////////////////////////////////////////////////
304
305 //---------------------------------------------------------
_Get_Difference(const CSG_Grid & Phi,int Level)306 bool CGridding_Spline_MBA_Grid::_Get_Difference(const CSG_Grid &Phi, int Level)
307 {
308 CSG_Simple_Statistics Differences;
309
310 double Scale = m_Points.Get_Cellsize() / Phi.Get_Cellsize();
311
312 for(int y=0; y<m_Points.Get_NY(); y++)
313 {
314 for(int x=0; x<m_Points.Get_NX(); x++)
315 {
316 if( !m_Points.is_NoData(x, y) )
317 {
318 double z = m_Points.asDouble(x, y) - BA_Get_Phi(Phi, x * Scale, y * Scale);
319
320 m_Points.Set_Value(x, y, z);
321
322 if( fabs(z) > m_Epsilon )
323 {
324 Differences += fabs(z);
325 }
326 }
327 }
328 }
329
330 //-----------------------------------------------------
331 Message_Fmt("\n%s:%d %s:%d %s:%f %s:%f",
332 _TL("level" ), Level + 1,
333 _TL("errors" ), (int)Differences.Get_Count (),
334 _TL("maximum"), Differences.Get_Maximum(),
335 _TL("mean" ), Differences.Get_Mean ()
336 );
337
338 Process_Set_Text(CSG_String::Format("%s %d [%d]", _TL("Level"), Level + 1, (int)Differences.Get_Count()));
339
340 return( Differences.Get_Maximum() > m_Epsilon );
341 }
342
343
344 ///////////////////////////////////////////////////////////
345 // //
346 ///////////////////////////////////////////////////////////
347
348 //---------------------------------------------------------
BA_Get_B(int i,double d) const349 inline double CGridding_Spline_MBA_Grid::BA_Get_B(int i, double d) const
350 {
351 switch( i )
352 {
353 case 0: d = 1. - d; return( d*d*d / 6. );
354
355 case 1: return( ( 3. * d*d*d - 6. * d*d + 4.) / 6. );
356
357 case 2: return( (-3. * d*d*d + 3. * d*d + 3. * d + 1.) / 6. );
358
359 case 3: return( d*d*d / 6. );
360
361 default: return( 0. );
362 }
363 }
364
365 //---------------------------------------------------------
BA_Set_Phi(CSG_Grid & Phi,double Cellsize)366 bool CGridding_Spline_MBA_Grid::BA_Set_Phi(CSG_Grid &Phi, double Cellsize)
367 {
368 int n = 4 + (int)(M_GET_MAX(m_pGrid->Get_XRange(), m_pGrid->Get_YRange()) / Cellsize);
369
370 Phi.Create(SG_DATATYPE_Float, n, n, Cellsize, m_pGrid->Get_XMin(), m_pGrid->Get_YMin());
371
372 CSG_Grid Delta(Phi.Get_System());
373
374 //-----------------------------------------------------
375 double Scale = m_Points.Get_Cellsize() / Phi.Get_Cellsize();
376
377 for(int yy=0; yy<m_Points.Get_NY(); yy++) for(int xx=0; xx<m_Points.Get_NX(); xx++)
378 {
379 if( m_Points.is_NoData(xx, yy) )
380 {
381 continue;
382 }
383
384 double p_x = xx * Scale; int x = (int)p_x;
385 double p_y = yy * Scale; int y = (int)p_y;
386 double p_z = m_Points.asDouble(xx, yy);
387
388 if( x >= 0 && x < Phi.Get_NX() - 3 && y >= 0 && y < Phi.Get_NY() - 3 )
389 {
390 int iy; double W[4][4], SW2 = 0.;
391
392 for(iy=0; iy<4; iy++) // compute W[k,l] and Sum[a=0-3, b=0-3](W�[a,b])
393 {
394 double wy = BA_Get_B(iy, p_y - y);
395
396 for(int ix=0; ix<4; ix++)
397 {
398 SW2 += SG_Get_Square(W[iy][ix] = wy * BA_Get_B(ix, p_x - x));
399 }
400 }
401
402 if( SW2 > 0. )
403 {
404 p_z /= SW2;
405
406 for(iy=0; iy<4; iy++)
407 {
408 for(int ix=0; ix<4; ix++)
409 {
410 double wxy = W[iy][ix];
411
412 Delta.Add_Value(x + ix, y + iy, wxy*wxy*wxy * p_z); // numerator
413 Phi .Add_Value(x + ix, y + iy, wxy*wxy ); // denominator
414 }
415 }
416 }
417 }
418 }
419
420 //-----------------------------------------------------
421 #pragma omp parallel for
422 for(int y=0; y<Phi.Get_NY(); y++)
423 {
424 for(int x=0; x<Phi.Get_NX(); x++)
425 {
426 double z = Phi.asDouble(x, y);
427
428 if( z != 0. )
429 {
430 Phi.Set_Value(x, y, Delta.asDouble(x, y) / z);
431 }
432 }
433 }
434
435 //-----------------------------------------------------
436 return( true );
437 }
438
439 //---------------------------------------------------------
BA_Get_Phi(const CSG_Grid & Phi,double px,double py) const440 double CGridding_Spline_MBA_Grid::BA_Get_Phi(const CSG_Grid &Phi, double px, double py) const
441 {
442 double z = 0.;
443
444 int x = (int)px;
445 int y = (int)py;
446
447 if( x >= 0 && x < Phi.Get_NX() - 3 && y >= 0 && y < Phi.Get_NY() - 3 )
448 {
449 for(int iy=0; iy<4; iy++)
450 {
451 double by = BA_Get_B(iy, py - y);
452
453 for(int ix=0; ix<4; ix++)
454 {
455 z += by * BA_Get_B(ix, px - x) * Phi.asDouble(x + ix, y + iy);
456 }
457 }
458 }
459
460 return( z );
461 }
462
463 //---------------------------------------------------------
BA_Set_Grid(const CSG_Grid & Phi,bool bAdd)464 void CGridding_Spline_MBA_Grid::BA_Set_Grid(const CSG_Grid &Phi, bool bAdd)
465 {
466 double d = m_pGrid->Get_Cellsize() / Phi.Get_Cellsize();
467
468 #pragma omp parallel for
469 for(int y=0; y<m_pGrid->Get_NY(); y++)
470 {
471 double py = d * y;
472
473 for(int x=0; x<m_pGrid->Get_NX(); x++)
474 {
475 double px = d * x;
476
477 if( bAdd )
478 { m_pGrid->Add_Value(x, y, BA_Get_Phi(Phi, px, py)); }
479 else
480 { m_pGrid->Set_Value(x, y, BA_Get_Phi(Phi, px, py)); }
481 }
482 }
483 }
484
485
486 ///////////////////////////////////////////////////////////
487 // //
488 // //
489 // //
490 ///////////////////////////////////////////////////////////
491
492 //---------------------------------------------------------
493