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_TPS_TIN.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_TPS_TIN.h"
52
53
54 ///////////////////////////////////////////////////////////
55 // //
56 // //
57 // //
58 ///////////////////////////////////////////////////////////
59
60 //---------------------------------------------------------
CGridding_Spline_TPS_TIN(void)61 CGridding_Spline_TPS_TIN::CGridding_Spline_TPS_TIN(void)
62 {
63 Set_Name (_TL("Thin Plate Spline (TIN)"));
64
65 Set_Author ("O.Conrad (c) 2006");
66
67 Set_Description (_TW(
68 "Creates a 'Thin Plate Spline' function for each triangle of a TIN "
69 "and uses it for subsequent gridding. The TIN is internally created "
70 "from the scattered data points input. The 'Neighbourhood' option "
71 "determines the number of points used for the spline generation. "
72 "'Immediate neighbourhood' includes the points of the triangle as well as the "
73 "immediate neighbour points. 'Level 1' adds the neighbours of the "
74 "immediate neighbourhood and 'level 2' adds the neighbours of 'level 1' "
75 "neighbours too. A higher neighbourhood degree reduces sharp breaks "
76 "but also increases the computation time. "
77 ));
78
79 Add_Reference(
80 "Donato G., Belongie S.", "2002",
81 "Approximation Methods for Thin Plate Spline Mappings and Principal Warps",
82 "In: Heyden, A., Sparr, G., Nielsen, M., Johansen, P. [Eds.]: "
83 "Computer Vision - ECCV 2002: 7th European Conference on Computer Vision, Copenhagen, Denmark, May 28-31, 2002, "
84 "Proceedings, Part III, Lecture Notes in Computer Science. Springer-Verlag Heidelberg; pp.21-31."
85 );
86
87 Add_Reference(
88 "Elonen, J.", "2005",
89 "Thin Plate Spline editor - an example program in C++",
90 "", SG_T("http://elonen.iki.fi/code/tpsdemo/index.html")
91 );
92
93 //-----------------------------------------------------
94 Parameters.Add_Double("",
95 "REGULARISATION" , _TL("Regularisation"),
96 _TL(""),
97 0.0001, 0., true
98 );
99
100 Parameters.Add_Choice("",
101 "LEVEL" , _TL("Neighbourhood"),
102 _TL(""),
103 CSG_String::Format("%s|%s|%s",
104 _TL("immediate"),
105 _TL("level 1"),
106 _TL("level 2")
107 ), 1
108 );
109
110 Parameters.Add_Bool("",
111 "FRAME" , _TL("Add Frame"),
112 _TL(""),
113 true
114 );
115 }
116
117
118 ///////////////////////////////////////////////////////////
119 // //
120 ///////////////////////////////////////////////////////////
121
122 //---------------------------------------------------------
_Initialise(void)123 bool CGridding_Spline_TPS_TIN::_Initialise(void)
124 {
125 m_Regularisation = Parameters("REGULARISATION")->asDouble();
126 m_Level = Parameters("LEVEL")->asInt();
127
128 m_Points = NULL;
129 m_nPoints_Buf = 0;
130
131 return( true );
132 }
133
134 //---------------------------------------------------------
_Finalise(void)135 bool CGridding_Spline_TPS_TIN::_Finalise(void)
136 {
137 if( m_Points )
138 {
139 SG_Free(m_Points);
140 }
141
142 m_Points = NULL;
143 m_nPoints_Buf = 0;
144
145 return( true );
146 }
147
148
149 ///////////////////////////////////////////////////////////
150 // //
151 ///////////////////////////////////////////////////////////
152
153 //---------------------------------------------------------
On_Execute(void)154 bool CGridding_Spline_TPS_TIN::On_Execute(void)
155 {
156 bool bResult = false;
157
158 CSG_TIN TIN;
159
160 if( Initialize() && _Initialise() && _Get_TIN(TIN) )
161 {
162 for(int i=0; i<TIN.Get_Triangle_Count() && Set_Progress(i, TIN.Get_Triangle_Count()); i++)
163 {
164 _Set_Triangle(TIN.Get_Triangle(i));
165 }
166
167 _Finalise();
168
169 bResult = true;
170 }
171
172 return( bResult );
173 }
174
175
176 ///////////////////////////////////////////////////////////
177 // //
178 ///////////////////////////////////////////////////////////
179
180 //---------------------------------------------------------
_Set_Triangle(CSG_TIN_Triangle * pTriangle)181 void CGridding_Spline_TPS_TIN::_Set_Triangle(CSG_TIN_Triangle *pTriangle)
182 {
183 if( m_pGrid->Get_Extent().Intersects(pTriangle->Get_Extent()) != INTERSECTION_None )
184 {
185 int i, j;
186
187 for(j=0, m_nPoints=0; j<3; j++)
188 {
189 CSG_TIN_Node *pPoint = pTriangle->Get_Node(j);
190
191 for(i=0; i<pPoint->Get_Neighbor_Count(); i++)
192 {
193 _Add_Points(pPoint->Get_Neighbor(i), 0);
194 }
195 }
196
197 CSG_Thin_Plate_Spline Spline;
198
199 for(i=0; i<m_nPoints; i++)
200 {
201 CSG_TIN_Node *pPoint = m_Points[i];
202
203 Spline.Add_Point(pPoint->Get_Point().x, pPoint->Get_Point().y, pPoint->asDouble(0));
204 }
205
206 if( Spline.Create(m_Regularisation, true) )
207 {
208 _Set_Grid(pTriangle, Spline);
209 }
210 }
211 }
212
213 //---------------------------------------------------------
_Set_Grid(CSG_TIN_Triangle * pTriangle,CSG_Thin_Plate_Spline & Spline)214 void CGridding_Spline_TPS_TIN::_Set_Grid(CSG_TIN_Triangle *pTriangle, CSG_Thin_Plate_Spline &Spline)
215 {
216 int ix, iy, ax, ay, bx, by;
217 double x, y, xMin, yMin;
218
219 ax = m_pGrid->Get_System().Get_xWorld_to_Grid(pTriangle->Get_Extent().Get_XMin()); if( ax < 0 ) ax = 0;
220 ay = m_pGrid->Get_System().Get_yWorld_to_Grid(pTriangle->Get_Extent().Get_YMin()); if( ay < 0 ) ay = 0;
221 bx = m_pGrid->Get_System().Get_xWorld_to_Grid(pTriangle->Get_Extent().Get_XMax()); if( bx >= m_pGrid->Get_NX() - 1 ) bx = m_pGrid->Get_NX() - 2;
222 by = m_pGrid->Get_System().Get_yWorld_to_Grid(pTriangle->Get_Extent().Get_YMax()); if( by >= m_pGrid->Get_NY() - 1 ) by = m_pGrid->Get_NY() - 2;
223 xMin = m_pGrid->Get_System().Get_xGrid_to_World(ax);
224 yMin = m_pGrid->Get_System().Get_yGrid_to_World(ay);
225
226 for(iy=ay, y=yMin; iy<=by; iy++, y+=m_pGrid->Get_Cellsize())
227 {
228 for(ix=ax, x=xMin; ix<=bx; ix++, x+=m_pGrid->Get_Cellsize())
229 {
230 if( pTriangle->is_Containing(x, y) )
231 {
232 m_pGrid->Set_Value(ix, iy, Spline.Get_Value(x, y));
233 }
234 }
235 }
236 }
237
238
239 ///////////////////////////////////////////////////////////
240 // //
241 ///////////////////////////////////////////////////////////
242
243 //---------------------------------------------------------
_Add_Points(CSG_TIN_Node * pPoint,int Level)244 void CGridding_Spline_TPS_TIN::_Add_Points(CSG_TIN_Node *pPoint, int Level)
245 {
246 _Add_Point(pPoint);
247
248 if( Level < m_Level )
249 {
250 for(int j=0; j<pPoint->Get_Neighbor_Count(); j++)
251 {
252 CSG_TIN_Node *p = pPoint->Get_Neighbor(j);
253
254 for(int i=0; i<p->Get_Neighbor_Count(); i++)
255 {
256 _Add_Points(p->Get_Neighbor(i), Level + 1);
257 }
258 }
259 }
260 }
261
262 //---------------------------------------------------------
_Add_Point(CSG_TIN_Node * pPoint)263 bool CGridding_Spline_TPS_TIN::_Add_Point(CSG_TIN_Node *pPoint)
264 {
265 for(int i=0; i<m_nPoints; i++)
266 {
267 if( m_Points[i] == pPoint )
268 {
269 return( false );
270 }
271 }
272
273 if( m_nPoints_Buf <= m_nPoints )
274 {
275 m_nPoints_Buf += 16;
276 m_Points = (CSG_TIN_Node **)SG_Realloc(m_Points, m_nPoints_Buf * sizeof(CSG_TIN_Node *));
277 }
278
279 m_Points[m_nPoints++] = pPoint;
280
281 return( true );
282 }
283
284
285 ///////////////////////////////////////////////////////////
286 // //
287 ///////////////////////////////////////////////////////////
288
289 //---------------------------------------------------------
_Get_TIN(CSG_TIN & TIN)290 bool CGridding_Spline_TPS_TIN::_Get_TIN(CSG_TIN &TIN)
291 {
292 TIN.Destroy();
293
294 bool bFrame = Parameters("FRAME" )->asBool ();
295 int zField = Parameters("FIELD" )->asInt ();
296 CSG_Shapes *pShapes = Parameters("SHAPES")->asShapes();
297
298 double x[4], y[4], z[4], dMin[4];
299
300 x[0] = m_pGrid->Get_Extent().Get_XMin(); y[0] = m_pGrid->Get_Extent().Get_YMin(); dMin[0] = -1.;
301 x[1] = m_pGrid->Get_Extent().Get_XMin(); y[1] = m_pGrid->Get_Extent().Get_YMax(); dMin[1] = -1.;
302 x[2] = m_pGrid->Get_Extent().Get_XMax(); y[2] = m_pGrid->Get_Extent().Get_YMax(); dMin[2] = -1.;
303 x[3] = m_pGrid->Get_Extent().Get_XMax(); y[3] = m_pGrid->Get_Extent().Get_YMin(); dMin[3] = -1.;
304
305 TIN.Add_Field("Z", pShapes->Get_Field_Type(zField));
306
307 for(int iShape=0; iShape<pShapes->Get_Count(); iShape++)
308 {
309 CSG_Shape *pShape = pShapes->Get_Shape(iShape);
310
311 if( !pShape->is_NoData(zField) )
312 {
313 for(int iPart=0; iPart<pShape->Get_Part_Count(); iPart++)
314 {
315 for(int iPoint=0; iPoint<pShape->Get_Point_Count(iPart); iPoint++)
316 {
317 TSG_Point p = pShape->Get_Point(iPoint, iPart);
318
319 TIN.Add_Node(p, NULL, false)->Set_Value(0, pShape->asDouble(zField));
320
321 if( bFrame )
322 {
323 for(int iCorner=0; iCorner<4; iCorner++)
324 {
325 double d = SG_Get_Distance(p.x, p.y, x[iCorner], y[iCorner]);
326
327 if( dMin[iCorner] < 0. || d < dMin[iCorner] )
328 {
329 dMin[iCorner] = d;
330 z [iCorner] = pShape->asDouble(zField);
331 }
332 }
333 }
334 }
335 }
336 }
337 }
338
339 //-----------------------------------------------------
340 if( bFrame )
341 {
342 for(int iCorner=0; iCorner<4; iCorner++)
343 {
344 if( dMin[iCorner] >= 0. )
345 {
346 CSG_Point p(x[iCorner], y[iCorner]);
347
348 TIN.Add_Node(p, NULL, false)->Set_Value(0, z[iCorner]);
349 }
350 }
351 }
352
353 TIN.Update();
354
355 return( TIN.Get_Triangle_Count() > 0 );
356 }
357
358
359 ///////////////////////////////////////////////////////////
360 // //
361 // //
362 // //
363 ///////////////////////////////////////////////////////////
364
365 //---------------------------------------------------------
366