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_3D.cpp //
14 // //
15 // Copyright (C) 2019 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 Hamburg //
44 // Germany //
45 // //
46 ///////////////////////////////////////////////////////////
47
48 //---------------------------------------------------------
49 #include "Gridding_Spline_MBA_3D.h"
50
51
52 ///////////////////////////////////////////////////////////
53 // //
54 // //
55 // //
56 ///////////////////////////////////////////////////////////
57
58 //---------------------------------------------------------
CGridding_Spline_MBA_3D(void)59 CGridding_Spline_MBA_3D::CGridding_Spline_MBA_3D(void)
60 {
61 Set_Name (_TL("Multilevel B-Spline (3D)"));
62
63 Set_Author ("O.Conrad (c) 2019");
64
65 Set_Description (_TW(
66 "Multilevel B-spline algorithm for spatial interpolation of scattered data "
67 "as proposed by Lee, Wolberg and Shin (1997) modified for 3D data.\n"
68 "The algorithm makes use of a coarse-to-fine hierarchy of control lattices to "
69 "generate a sequence of bicubic B-spline functions, whose sum approaches the "
70 "desired interpolation function. Performance gains are realized by using "
71 "B-spline refinement to reduce the sum of these functions into one equivalent "
72 "B-spline function. "
73 "\n\n"
74 "The 'Maximum Level' determines the maximum size of the final B-spline matrix "
75 "and increases exponential with each level. Where level=10 requires about 1mb "
76 "level=12 needs about 16mb and level=14 about 256mb(!) of additional memory. "
77 ));
78
79 Add_Reference(
80 "Lee, S., Wolberg, G., Shin, S.Y.", "1997",
81 "Scattered Data Interpolation with Multilevel B-Splines",
82 "IEEE Transactions On Visualisation And Computer Graphics, Vol.3, No.3., p.228-244.",
83 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"),
84 SG_T("ResearchGate")
85 );
86
87 //-----------------------------------------------------
88 Parameters.Add_Shapes("",
89 "POINTS" , _TL("Points"),
90 _TL(""),
91 PARAMETER_INPUT, SHAPE_TYPE_Point
92 );
93
94 Parameters.Add_Table_Field("POINTS",
95 "Z_FIELD" , _TL("Z"),
96 _TL("")
97 );
98
99 Parameters.Add_Double("POINTS",
100 "Z_SCALE" , _TL("Z Factor"),
101 _TL(""),
102 1., 0., true
103 );
104
105 Parameters.Add_Table_Field("POINTS",
106 "V_FIELD" , _TL("Value"),
107 _TL("")
108 );
109
110 //-----------------------------------------------------
111 m_Grid_Target.Create(&Parameters, false, "", "TARGET_");
112
113 m_Grid_Target.Add_Grids("GRIDS", _TL("Grid Collection"), false, true);
114
115 //-----------------------------------------------------
116 Parameters.Add_Double(
117 "", "EPSILON" , _TL("Threshold Error"),
118 _TL(""),
119 0.0001, 0., true
120 );
121
122 Parameters.Add_Int(
123 "", "LEVEL_MAX" , _TL("Maximum Level"),
124 _TL(""),
125 11, 1, true, 14, true
126 );
127 }
128
129
130 ///////////////////////////////////////////////////////////
131 // //
132 ///////////////////////////////////////////////////////////
133
134 //---------------------------------------------------------
On_Parameter_Changed(CSG_Parameters * pParameters,CSG_Parameter * pParameter)135 int CGridding_Spline_MBA_3D::On_Parameter_Changed(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
136 {
137 if( pParameter->Cmp_Identifier("POINTS") )
138 {
139 m_Grid_Target.Set_User_Defined(pParameters, pParameter->asShapes());
140 }
141
142 if( pParameter->Cmp_Identifier("POINTS") || pParameter->Cmp_Identifier("Z_FIELD") )
143 {
144 CSG_Shapes *pPoints = (*pParameters)("POINTS")->asShapes();
145
146 if( pPoints )
147 {
148 int zField = pPoints->Get_Vertex_Type() == SG_VERTEX_TYPE_XY ? (*pParameters)("Z_FIELD")->asInt() : -1;
149
150 m_Grid_Target.Set_User_Defined_ZLevels(pParameters,
151 zField < 0 ? pPoints->Get_ZMin() : pPoints->Get_Minimum(zField),
152 zField < 0 ? pPoints->Get_ZMax() : pPoints->Get_Maximum(zField), 10
153 );
154 }
155 }
156
157 m_Grid_Target.On_Parameter_Changed(pParameters, pParameter);
158
159 return( CSG_Tool::On_Parameter_Changed(pParameters, pParameter) );
160 }
161
162 //---------------------------------------------------------
On_Parameters_Enable(CSG_Parameters * pParameters,CSG_Parameter * pParameter)163 int CGridding_Spline_MBA_3D::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
164 {
165 if( pParameter->Cmp_Identifier("POINTS") )
166 {
167 pParameters->Set_Enabled("Z_FIELD", pParameter->asShapes() && pParameter->asShapes()->Get_Vertex_Type() == SG_VERTEX_TYPE_XY);
168 }
169
170 m_Grid_Target.On_Parameters_Enable(pParameters, pParameter);
171
172 return( CSG_Tool::On_Parameters_Enable(pParameters, pParameter) );
173 }
174
175
176 ///////////////////////////////////////////////////////////
177 // //
178 ///////////////////////////////////////////////////////////
179
180 //---------------------------------------------------------
On_Execute(void)181 bool CGridding_Spline_MBA_3D::On_Execute(void)
182 {
183 bool bResult = false;
184
185 if( !Initialize() )
186 {
187 return( false );
188 }
189
190 m_Epsilon = Parameters("EPSILON")->asDouble();
191
192 double Cellsize = M_GET_MAX(M_GET_MAX(m_pGrids->Get_XRange(), m_pGrids->Get_YRange()), m_pGrids->Get_ZRange());
193
194 bResult = _Set_MBA(Cellsize);
195
196 m_Points.Destroy();
197
198 if( m_zField >= 0 && m_zField != m_pGrids->Get_Z_Attribute() ) // z factor != 1 !!
199 {
200 int zField = m_pGrids->Get_Z_Attribute();
201
202 m_pGrids->Set_Z_Attribute (m_zField);
203 m_pGrids->Set_Z_Name_Field(m_zField);
204
205 m_pGrids->Del_Attribute(zField);
206 }
207
208 Finalize();
209
210 return( bResult );
211 }
212
213
214 ///////////////////////////////////////////////////////////
215 // //
216 ///////////////////////////////////////////////////////////
217
218 //---------------------------------------------------------
Initialize(void)219 bool CGridding_Spline_MBA_3D::Initialize(void)
220 {
221 CSG_Shapes *pPoints = Parameters("POINTS")->asShapes();
222
223 int zField = pPoints->Get_Vertex_Type() == SG_VERTEX_TYPE_XY ? Parameters("Z_FIELD")->asInt() : -1;
224
225 int vField = Parameters("V_FIELD")->asInt();
226
227 if( (m_pGrids = m_Grid_Target.Get_Grids("GRIDS")) == NULL )
228 {
229 return( false );
230 }
231
232 m_pGrids->Fmt_Name("%s.%s [%s]", pPoints->Get_Name(), Parameters("V_FIELD")->asString(), Get_Name().c_str());
233
234 m_zCellsize = Parameters("TARGET_" "USER_ZSIZE")->asDouble();
235
236 double zScale = Parameters("Z_SCALE")->asDouble();
237
238 if( zScale == 0. )
239 {
240 Error_Set(_TL("Z factor is zero! Please use 2D instead of 3D interpolation."));
241
242 return( false );
243 }
244
245 m_zField = zScale == 1. ? -1 : m_pGrids->Get_Z_Attribute();
246
247 if( m_zField >= 0 )
248 {
249 m_pGrids->Add_Attribute("Z_Scaled", SG_DATATYPE_Double);
250
251 for(int iz=0; iz<m_pGrids->Get_NZ(); iz++)
252 {
253 m_pGrids->Get_Attributes(iz).Set_Value("Z_Scaled", m_pGrids->Get_Z(iz) * zScale);
254 }
255
256 m_pGrids->Set_Z_Attribute(m_pGrids->Get_Attributes().Get_Field_Count() - 1);
257
258 m_zCellsize *= zScale;
259 }
260
261 //-----------------------------------------------------
262 m_Points.Destroy();
263
264 for(int i=0; i<pPoints->Get_Count(); i++)
265 {
266 CSG_Shape *pPoint = pPoints->Get_Shape(i);
267
268 if( (zField < 0 || !pPoint->is_NoData(zField)) && !pPoint->is_NoData(vField) )
269 {
270 CSG_Vector p(4);
271
272 p[0] = pPoint->Get_Point(0).x;
273 p[1] = pPoint->Get_Point(0).y;
274 p[2] = zScale * (zField < 0 ? pPoint->Get_Z(0)
275 : pPoint->asDouble(zField));
276 p[3] = pPoint->asDouble(vField) - pPoints->Get_Mean(vField); // detrend!
277
278 m_Points.Add_Row(p);
279 }
280 }
281
282 return( m_Points.Get_NRows() > 0 );
283 }
284
285 //---------------------------------------------------------
Finalize(void)286 bool CGridding_Spline_MBA_3D::Finalize(void)
287 {
288 double Mean = Parameters("POINTS")->asShapes()->Get_Mean(Parameters("V_FIELD")->asInt());
289
290 if( Mean ) // detrend!
291 {
292 m_pGrids->Add(Mean);
293 }
294
295 return( true );
296 }
297
298
299 ///////////////////////////////////////////////////////////
300 // //
301 ///////////////////////////////////////////////////////////
302
303 //---------------------------------------------------------
_Set_MBA(double Cellsize)304 bool CGridding_Spline_MBA_3D::_Set_MBA(double Cellsize)
305 {
306 CSG_Grids Phi;
307
308 bool bContinue = true;
309
310 int Levels = Parameters("LEVEL_MAX")->asInt();
311
312 for(int Level=0; bContinue && Level<Levels && Process_Get_Okay(false); Level++, Cellsize/=2.)
313 {
314 bContinue = BA_Set_Phi(Phi, Cellsize) && _Get_Difference(Phi, Level);
315
316 BA_Set_Grids(Phi, Level > 0);
317 }
318
319 return( true );
320 }
321
322
323 ///////////////////////////////////////////////////////////
324 // //
325 ///////////////////////////////////////////////////////////
326
327 //---------------------------------------------------------
_Get_Difference(const CSG_Grids & Phi,int Level)328 bool CGridding_Spline_MBA_3D::_Get_Difference(const CSG_Grids &Phi, int Level)
329 {
330 CSG_Simple_Statistics Differences;
331
332 for(int i=0; i<m_Points.Get_NRows(); i++)
333 {
334 CSG_Vector p(4, m_Points[i]);
335
336 p[0] = (p[0] - Phi.Get_XMin()) / Phi.Get_Cellsize();
337 p[1] = (p[1] - Phi.Get_YMin()) / Phi.Get_Cellsize();
338 p[2] = (p[2] - Phi.Get_ZMin()) / Phi.Get_Cellsize();
339 p[3] = p[3] - BA_Get_Phi(Phi, p[0], p[1], p[2]);
340
341 m_Points[i][3] = p[3];
342
343 if( fabs(p[3]) > m_Epsilon )
344 {
345 Differences += fabs(p[3]);
346 }
347 }
348
349 //-----------------------------------------------------
350 Message_Fmt("\n%s:%d %s:%d %s:%f %s:%f",
351 _TL("level" ), Level + 1,
352 _TL("errors" ), (int)Differences.Get_Count (),
353 _TL("maximum"), Differences.Get_Maximum(),
354 _TL("mean" ), Differences.Get_Mean ()
355 );
356
357 Process_Set_Text(CSG_String::Format("%s %d [%d]", _TL("Level"), Level + 1, (int)Differences.Get_Count()));
358
359 return( Differences.Get_Maximum() > m_Epsilon );
360 }
361
362
363 ///////////////////////////////////////////////////////////
364 // //
365 ///////////////////////////////////////////////////////////
366
367 //---------------------------------------------------------
BA_Get_B(int i,double d) const368 inline double CGridding_Spline_MBA_3D::BA_Get_B(int i, double d) const
369 {
370 switch( i )
371 {
372 case 0: d = 1. - d; return( d*d*d / 6. );
373
374 case 1: return( ( 3. * d*d*d - 6. * d*d + 4.) / 6. );
375
376 case 2: return( (-3. * d*d*d + 3. * d*d + 3. * d + 1.) / 6. );
377
378 case 3: return( d*d*d / 6. );
379
380 default: return( 0. );
381 }
382 }
383
384 //---------------------------------------------------------
BA_Set_Phi(CSG_Grids & Phi,double Cellsize)385 bool CGridding_Spline_MBA_3D::BA_Set_Phi(CSG_Grids &Phi, double Cellsize)
386 {
387 int n = 4 + (int)(M_GET_MAX(M_GET_MAX(m_pGrids->Get_XRange(), m_pGrids->Get_YRange()), m_pGrids->Get_ZRange()) / Cellsize);
388
389 Phi.Create(n, n, n, Cellsize, m_pGrids->Get_XMin(), m_pGrids->Get_YMin(), m_pGrids->Get_ZMin(), SG_DATATYPE_Float);
390
391 CSG_Grids Delta(n, n, n, Cellsize, m_pGrids->Get_XMin(), m_pGrids->Get_YMin(), m_pGrids->Get_ZMin(), SG_DATATYPE_Float);
392
393 if( Phi.Get_NZ() < n || Delta.Get_NZ() < n )
394 {
395 Message_Fmt("\n%s", _TL("failed to allocate memory for phi calculation"));
396
397 return( false );
398 }
399
400 //-----------------------------------------------------
401 for(int i=0; i<m_Points.Get_NRows(); i++)
402 {
403 CSG_Vector p(4, m_Points[i]);
404
405 int x = (int)(p[0] = (p[0] - Phi.Get_XMin()) / Phi.Get_Cellsize());
406 int y = (int)(p[1] = (p[1] - Phi.Get_YMin()) / Phi.Get_Cellsize());
407 int z = (int)(p[2] = (p[2] - Phi.Get_ZMin()) / Phi.Get_Cellsize());
408
409 if( x >= 0 && x < Phi.Get_NX() - 3 && y >= 0 && y < Phi.Get_NY() - 3 && z >= 0 && z < Phi.Get_NZ() - 3 )
410 {
411 int iz; double W[4][4][4], SW2 = 0.;
412
413 for(iz=0; iz<4; iz++) // compute W[k,l] and Sum[a=0-3, b=0-3](W�[a,b])
414 {
415 double wz = BA_Get_B(iz, p[2] - z);
416
417 for(int iy=0; iy<4; iy++)
418 {
419 double wyz = wz * BA_Get_B(iy, p[1] - y);
420
421 for(int ix=0; ix<4; ix++)
422 {
423 SW2 += SG_Get_Square(W[iz][iy][ix] = wyz * BA_Get_B(ix, p[0] - x));
424 }
425 }
426 }
427
428 if( SW2 > 0. )
429 {
430 double dz = p[3] / SW2;
431
432 for(iz=0; iz<4; iz++)
433 {
434 for(int iy=0; iy<4; iy++)
435 {
436 for(int ix=0; ix<4; ix++)
437 {
438 double wxyz = W[iz][iy][ix];
439
440 Delta.Add_Value(x + ix, y + iy, z + iz, wxyz*wxyz*wxyz * dz); // numerator
441 Phi .Add_Value(x + ix, y + iy, z + iz, wxyz*wxyz ); // denominator
442 }
443 }
444 }
445 }
446 }
447 }
448
449 //-----------------------------------------------------
450 #pragma omp parallel for
451 for(int z=0; z<Phi.Get_NZ(); z++)
452 {
453 for(int y=0; y<Phi.Get_NY(); y++)
454 {
455 for(int x=0; x<Phi.Get_NX(); x++)
456 {
457 double v = Phi.asDouble(x, y, z);
458
459 if( v != 0. )
460 {
461 Phi.Set_Value(x, y, z, Delta.asDouble(x, y, z) / v);
462 }
463 }
464 }
465 }
466
467 //-----------------------------------------------------
468 return( true );
469 }
470
471 //---------------------------------------------------------
BA_Get_Phi(const CSG_Grids & Phi,double px,double py,double pz) const472 double CGridding_Spline_MBA_3D::BA_Get_Phi(const CSG_Grids &Phi, double px, double py, double pz) const
473 {
474 double v = 0.;
475
476 int x = (int)px; px -= x;
477 int y = (int)py; py -= y;
478 int z = (int)pz; pz -= z;
479
480 if( x >= 0 && x < Phi.Get_NX() - 3 && y >= 0 && y < Phi.Get_NY() - 3 && z >= 0 && z < Phi.Get_NZ() - 3 )
481 {
482 for(int iz=0; iz<4; iz++)
483 {
484 double bz = BA_Get_B(iz, pz);
485
486 for(int iy=0; iy<4; iy++)
487 {
488 double byz = bz * BA_Get_B(iy, py);
489
490 for(int ix=0; ix<4; ix++)
491 {
492 v += byz * BA_Get_B(ix, px) * Phi.asDouble(x + ix, y + iy, z + iz);
493 }
494 }
495 }
496 }
497
498 return( v );
499 }
500
501 //---------------------------------------------------------
BA_Set_Grids(const CSG_Grids & Phi,bool bAdd)502 void CGridding_Spline_MBA_3D::BA_Set_Grids(const CSG_Grids &Phi, bool bAdd)
503 {
504 double d = m_pGrids->Get_Cellsize() / Phi.Get_Cellsize();
505
506 #pragma omp parallel for
507 for(int z=0; z<m_pGrids->Get_NZ(); z++)
508 {
509 double pz = z * m_zCellsize / Phi.Get_Cellsize();
510
511 for(int y=0; y<m_pGrids->Get_NY(); y++)
512 {
513 double py = y * d;
514
515 for(int x=0; x<m_pGrids->Get_NX(); x++)
516 {
517 double px = x * d;
518
519 if( bAdd )
520 { m_pGrids->Add_Value(x, y, z, BA_Get_Phi(Phi, px, py, pz)); }
521 else
522 { m_pGrids->Set_Value(x, y, z, BA_Get_Phi(Phi, px, py, pz)); }
523 }
524 }
525 }
526 }
527
528
529 ///////////////////////////////////////////////////////////
530 // //
531 // //
532 // //
533 ///////////////////////////////////////////////////////////
534
535 //---------------------------------------------------------
536