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_BA.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_BA.h"
52
53
54 ///////////////////////////////////////////////////////////
55 // //
56 // //
57 // //
58 ///////////////////////////////////////////////////////////
59
60 //---------------------------------------------------------
CGridding_Spline_BA(void)61 CGridding_Spline_BA::CGridding_Spline_BA(void)
62 : CGridding_Spline_Base()
63 {
64 Set_Name (_TL("B-Spline Approximation"));
65
66 Set_Author ("O.Conrad (c) 2006");
67
68 Set_Description (_TW(
69 "Calculates B-spline functions for chosen level of detail. "
70 "This tool serves as the basis for the 'Multilevel B-spline Interpolation' "
71 "and is not suited as is for spatial data interpolation from "
72 "scattered data. "
73 ));
74
75 Add_Reference(
76 "Lee, S., Wolberg, G., Shin, S.Y.", "1997",
77 "Scattered Data Interpolation with Multilevel B-Splines",
78 "IEEE Transactions On Visualisation And Computer Graphics, Vol.3, No.3., p.228-244.",
79 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"),
80 SG_T("ResearchGate")
81 );
82
83 //-----------------------------------------------------
84 Parameters.Add_Double(
85 "", "LEVEL" , _TL("Range"),
86 _TL("B-spline range expressed as number of cells."),
87 1, 0.001, true
88 );
89 }
90
91
92 ///////////////////////////////////////////////////////////
93 // //
94 ///////////////////////////////////////////////////////////
95
96 //---------------------------------------------------------
On_Execute(void)97 bool CGridding_Spline_BA::On_Execute(void)
98 {
99 bool bResult = false;
100
101 if( Initialize(m_Points, true) )
102 {
103 double Cellsize = m_pGrid->Get_Cellsize() * Parameters("LEVEL")->asDouble();
104
105 CSG_Grid Phi;
106
107 if( BA_Set_Phi(Phi, Cellsize) )
108 {
109 BA_Set_Grid(Phi);
110
111 bResult = true;
112 }
113 }
114
115 m_Points.Clear();
116
117 return( bResult );
118 }
119
120
121 ///////////////////////////////////////////////////////////
122 // //
123 ///////////////////////////////////////////////////////////
124
125 //---------------------------------------------------------
BA_Get_B(int i,double d) const126 inline double CGridding_Spline_BA::BA_Get_B(int i, double d) const
127 {
128 switch( i )
129 {
130 case 0: d = 1. - d; return( d*d*d / 6. );
131
132 case 1: return( ( 3. * d*d*d - 6. * d*d + 4.) / 6. );
133
134 case 2: return( (-3. * d*d*d + 3. * d*d + 3. * d + 1.) / 6. );
135
136 case 3: return( d*d*d / 6. );
137
138 default: return( 0. );
139 }
140 }
141
142 //---------------------------------------------------------
BA_Set_Phi(CSG_Grid & Phi,double Cellsize)143 bool CGridding_Spline_BA::BA_Set_Phi(CSG_Grid &Phi, double Cellsize)
144 {
145 int nx = (int)((m_pGrid->Get_XRange()) / Cellsize);
146 int ny = (int)((m_pGrid->Get_YRange()) / Cellsize);
147
148 Phi.Create(SG_DATATYPE_Float, nx + 4, ny + 4, Cellsize, m_pGrid->Get_XMin(), m_pGrid->Get_YMin());
149
150 CSG_Grid Delta(Phi.Get_System());
151
152 //-----------------------------------------------------
153 for(int i=0; i<m_Points.Get_Count(); i++)
154 {
155 TSG_Point_Z p = m_Points[i];
156
157 int x = (int)(p.x = (p.x - Phi.Get_XMin()) / Phi.Get_Cellsize());
158 int y = (int)(p.y = (p.y - Phi.Get_YMin()) / Phi.Get_Cellsize());
159
160 if( x >= 0 && x < Phi.Get_NX() - 3 && y >= 0 && y < Phi.Get_NY() - 3 )
161 {
162 int iy; double W[4][4], SW2 = 0.;
163
164 for(iy=0; iy<4; iy++) // compute W[k,l] and Sum[a=0-3, b=0-3](W�[a,b])
165 {
166 double wy = BA_Get_B(iy, p.y - y);
167
168 for(int ix=0; ix<4; ix++)
169 {
170 SW2 += SG_Get_Square(W[iy][ix] = wy * BA_Get_B(ix, p.x - x));
171 }
172 }
173
174 if( SW2 > 0. )
175 {
176 p.z /= SW2;
177
178 for(iy=0; iy<4; iy++)
179 {
180 for(int ix=0; ix<4; ix++)
181 {
182 double wxy = W[iy][ix];
183
184 Delta.Add_Value(x + ix, y + iy, wxy*wxy*wxy * p.z); // numerator
185 Phi .Add_Value(x + ix, y + iy, wxy*wxy ); // denominator
186 }
187 }
188 }
189 }
190 }
191
192 //-----------------------------------------------------
193 #pragma omp parallel for
194 for(int y=0; y<Phi.Get_NY(); y++)
195 {
196 for(int x=0; x<Phi.Get_NX(); x++)
197 {
198 double z = Phi.asDouble(x, y);
199
200 if( z != 0. )
201 {
202 Phi.Set_Value(x, y, Delta.asDouble(x, y) / z);
203 }
204 }
205 }
206
207 //-----------------------------------------------------
208 return( true );
209 }
210
211
212 ///////////////////////////////////////////////////////////
213 // //
214 ///////////////////////////////////////////////////////////
215
216 //---------------------------------------------------------
BA_Get_Phi(const CSG_Grid & Phi,double px,double py) const217 double CGridding_Spline_BA::BA_Get_Phi(const CSG_Grid &Phi, double px, double py) const
218 {
219 double z = 0.;
220
221 int x = (int)px;
222 int y = (int)py;
223
224 if( x >= 0 && x < Phi.Get_NX() - 3 && y >= 0 && y < Phi.Get_NY() - 3 )
225 {
226 for(int iy=0; iy<4; iy++)
227 {
228 double by = BA_Get_B(iy, py - y);
229
230 for(int ix=0; ix<4; ix++)
231 {
232 z += by * BA_Get_B(ix, px - x) * Phi.asDouble(x + ix, y + iy);
233 }
234 }
235 }
236
237 return( z );
238 }
239
240 //---------------------------------------------------------
BA_Set_Grid(const CSG_Grid & Phi,bool bAdd)241 void CGridding_Spline_BA::BA_Set_Grid(const CSG_Grid &Phi, bool bAdd)
242 {
243 double d = m_pGrid->Get_Cellsize() / Phi.Get_Cellsize();
244
245 #pragma omp parallel for
246 for(int y=0; y<m_pGrid->Get_NY(); y++)
247 {
248 double py = d * y;
249
250 for(int x=0; x<m_pGrid->Get_NX(); x++)
251 {
252 double px = d * x;
253
254 if( bAdd )
255 { m_pGrid->Add_Value(x, y, BA_Get_Phi(Phi, px, py)); }
256 else
257 { m_pGrid->Set_Value(x, y, BA_Get_Phi(Phi, px, py)); }
258 }
259 }
260 }
261
262
263 ///////////////////////////////////////////////////////////
264 // //
265 // //
266 // //
267 ///////////////////////////////////////////////////////////
268
269 //---------------------------------------------------------
270