1 | /* |
---|
2 | SWWReader |
---|
3 | |
---|
4 | Reader of SWW files from within OpenSceneGraph. |
---|
5 | |
---|
6 | An SWW file is the visualization output of the pyVolution |
---|
7 | Shallow Water Wave equation solver. SWW files contain bedslope |
---|
8 | geometry and a sequence of stages at known timesteps. The |
---|
9 | internal format is NetCDF: |
---|
10 | |
---|
11 | netcdf demo.sww { |
---|
12 | dimensions: |
---|
13 | number_of_volumes = <integer> ; |
---|
14 | number_of_vertices = 3 ; |
---|
15 | number_of_points = <integer> ; |
---|
16 | number_of_timesteps = UNLIMITED ; |
---|
17 | variables: |
---|
18 | float x(number_of_points) ; |
---|
19 | float y(number_of_points) ; |
---|
20 | float z(number_of_points) ; |
---|
21 | int volumes(number_of_volumes, number_of_vertices) ; |
---|
22 | float time(number_of_timesteps) ; |
---|
23 | float stage(number_of_timesteps, number_of_points) ; |
---|
24 | |
---|
25 | copyright (C) 2004 Geoscience Australia |
---|
26 | */ |
---|
27 | |
---|
28 | #include <SWWReader.h> |
---|
29 | #include <string> |
---|
30 | #include <fstream> |
---|
31 | #include <netcdf.h> |
---|
32 | #include <osg/Notify> |
---|
33 | |
---|
34 | #include <stdlib.h> |
---|
35 | #include <stdio.h> |
---|
36 | |
---|
37 | #define DEF_CULLNEARZERO 0.001 |
---|
38 | #define DEF_CULLSTEEPANGLE 0.9 |
---|
39 | #define DEF_CULLSTART CULLNONE |
---|
40 | |
---|
41 | |
---|
42 | // only constructor, requires netcdf file |
---|
43 | SWWReader::SWWReader(const std::string& filename) |
---|
44 | { |
---|
45 | // assume failure until end of constructor |
---|
46 | _valid = false; |
---|
47 | |
---|
48 | // netcdf filename |
---|
49 | _filename = new std::string(filename); |
---|
50 | |
---|
51 | // netcdf open |
---|
52 | _status.push_back( nc_open(_filename->c_str(), NC_NOWRITE, &_ncid) ); |
---|
53 | if (this->_statusHasError()) return; |
---|
54 | |
---|
55 | // dimension ids |
---|
56 | _status.push_back( nc_inq_dimid(_ncid, "number_of_volumes", &_nvolumesid) ); |
---|
57 | _status.push_back( nc_inq_dimid(_ncid, "number_of_vertices", &_nverticesid) ); |
---|
58 | _status.push_back( nc_inq_dimid(_ncid, "number_of_points", &_npointsid) ); |
---|
59 | _status.push_back( nc_inq_dimid(_ncid, "number_of_timesteps", &_ntimestepsid) ); |
---|
60 | if (this->_statusHasError()) return; |
---|
61 | |
---|
62 | // dimension values |
---|
63 | _status.push_back( nc_inq_dimlen(_ncid, _nvolumesid, &_nvolumes) ); |
---|
64 | _status.push_back( nc_inq_dimlen(_ncid, _nverticesid, &_nvertices) ); |
---|
65 | _status.push_back( nc_inq_dimlen(_ncid, _npointsid, &_npoints) ); |
---|
66 | _status.push_back( nc_inq_dimlen(_ncid, _ntimestepsid, &_ntimesteps) ); |
---|
67 | if (this->_statusHasError()) return; |
---|
68 | |
---|
69 | // variable ids |
---|
70 | _status.push_back( nc_inq_varid (_ncid, "x", &_xid) ); |
---|
71 | _status.push_back( nc_inq_varid (_ncid, "y", &_yid) ); |
---|
72 | _status.push_back( nc_inq_varid (_ncid, "z", &_zid) ); |
---|
73 | _status.push_back( nc_inq_varid (_ncid, "volumes", &_volumesid) ); |
---|
74 | _status.push_back( nc_inq_varid (_ncid, "time", &_timeid) ); |
---|
75 | _status.push_back( nc_inq_varid (_ncid, "stage", &_stageid) ); |
---|
76 | if (this->_statusHasError()) return; |
---|
77 | |
---|
78 | // allocation of variable arrays, destructor responsible for cleanup |
---|
79 | _px = new float[_npoints]; |
---|
80 | _py = new float[_npoints]; |
---|
81 | _pz = new float[_npoints]; |
---|
82 | _ptime = new float[_ntimesteps]; |
---|
83 | _pvolumes = new unsigned int[_nvertices * _nvolumes]; |
---|
84 | _pstage = new float[_npoints]; |
---|
85 | |
---|
86 | // loading variables from netcdf file |
---|
87 | _status.push_back( nc_get_var_float (_ncid, _xid, _px) ); // x vertices |
---|
88 | _status.push_back( nc_get_var_float (_ncid, _yid, _py) ); // y vertices |
---|
89 | _status.push_back( nc_get_var_float (_ncid, _zid, _pz) ); // bedslope heights |
---|
90 | _status.push_back( nc_get_var_float (_ncid, _timeid, _ptime) ); // time array |
---|
91 | _status.push_back( nc_get_var_int (_ncid, _volumesid, (int *) _pvolumes) ); // triangle indices |
---|
92 | if (this->_statusHasError()) return; |
---|
93 | |
---|
94 | // loop index |
---|
95 | size_t iv; |
---|
96 | |
---|
97 | // vertex indices |
---|
98 | unsigned int v1index, v2index, v3index; |
---|
99 | |
---|
100 | // bedslope extents and resultant scale factors |
---|
101 | float xmin, xmax, xrange; |
---|
102 | float ymin, ymax, yrange; |
---|
103 | float zmin, zmax, zrange; |
---|
104 | xmin = _px[0]; |
---|
105 | xmax = _px[0]; |
---|
106 | ymin = _py[0]; |
---|
107 | ymax = _py[0]; |
---|
108 | zmin = _pz[0]; |
---|
109 | zmax = _pz[0]; |
---|
110 | for( iv=1; iv < _npoints; iv++ ) |
---|
111 | { |
---|
112 | if( _px[iv] < xmin ) xmin = _px[iv]; |
---|
113 | if( _px[iv] > xmax ) xmax = _px[iv]; |
---|
114 | if( _py[iv] < ymin ) ymin = _py[iv]; |
---|
115 | if( _py[iv] > ymax ) ymax = _py[iv]; |
---|
116 | if( _pz[iv] < zmin ) zmin = _pz[iv]; |
---|
117 | if( _pz[iv] > zmax ) zmax = _pz[iv]; |
---|
118 | } |
---|
119 | _xscale = 1.0/(xmax - xmin); _xoffset = xmin; |
---|
120 | _yscale = 1.0/(ymax - ymin); _yoffset = ymin; |
---|
121 | _zscale = (zmax==zmin) ? 1.0 : 1.0/(zmax - zmin); _zoffset = zmin; |
---|
122 | |
---|
123 | #ifdef DEBUG |
---|
124 | std::cout << "xmin: " << xmin << std::endl; |
---|
125 | std::cout << "xmax: " << xmax << std::endl; |
---|
126 | std::cout << "ymin: " << ymin << std::endl; |
---|
127 | std::cout << "ymax: " << ymax << std::endl; |
---|
128 | std::cout << "zmin: " << zmin << std::endl; |
---|
129 | std::cout << "zmax: " << zmax << std::endl; |
---|
130 | std::cout << "xscale: " << _xscale << std::endl; |
---|
131 | std::cout << "yscale: " << _yscale << std::endl; |
---|
132 | std::cout << "zscale: " << _zscale << std::endl; |
---|
133 | #endif |
---|
134 | |
---|
135 | // culling defaults |
---|
136 | _cullsetting = DEF_CULLSTART; |
---|
137 | _cullnearzero = DEF_CULLNEARZERO; |
---|
138 | _cullsteepangle = DEF_CULLSTEEPANGLE; |
---|
139 | |
---|
140 | // compute triangle connectivity, a list (indexed by vertex number) |
---|
141 | // of lists (indices of triangles sharing this vertex) |
---|
142 | _connectivity = std::vector<triangle_list>(_npoints); |
---|
143 | for (iv=0; iv < _nvolumes; iv++) |
---|
144 | { |
---|
145 | v1index = _pvolumes[3*iv+0]; |
---|
146 | v2index = _pvolumes[3*iv+1]; |
---|
147 | v3index = _pvolumes[3*iv+2]; |
---|
148 | _connectivity.at(v1index).push_back(iv); |
---|
149 | _connectivity.at(v2index).push_back(iv); |
---|
150 | _connectivity.at(v3index).push_back(iv); |
---|
151 | } |
---|
152 | |
---|
153 | // load bedslope vertex array, shifting and scaling vertices to unit cube |
---|
154 | // centred about the origin |
---|
155 | _bedslopevertices = new osg::Vec3Array; |
---|
156 | for (iv=0; iv < _npoints; iv++) |
---|
157 | _bedslopevertices->push_back( osg::Vec3( (_px[iv]-_xoffset)*_xscale - 0.5, |
---|
158 | (_py[iv]-_yoffset)*_yscale - 0.5, |
---|
159 | (_pz[iv]-_zoffset)*_zscale - 0.5) ); |
---|
160 | |
---|
161 | // load bedslope index array, pvolumes array indexes into x, y and z |
---|
162 | _bedslopeindices = new osg::UIntArray; |
---|
163 | for (iv=0; iv < _nvolumes*_nvertices; iv++) |
---|
164 | _bedslopeindices->push_back( _pvolumes[iv] ); |
---|
165 | |
---|
166 | // calculate bedslope primitive normal and centroid arrays |
---|
167 | _bedslopenormals = new osg::Vec3Array; |
---|
168 | _bedslopecentroids = new osg::Vec3Array; |
---|
169 | osg::Vec3 v1, v2, v3, side1, side2, nrm; |
---|
170 | for (iv=0; iv < _nvolumes; iv++) |
---|
171 | { |
---|
172 | v1index = _pvolumes[3*iv+0]; |
---|
173 | v2index = _pvolumes[3*iv+1]; |
---|
174 | v3index = _pvolumes[3*iv+2]; |
---|
175 | |
---|
176 | v1 = _bedslopevertices->at(v1index); |
---|
177 | v2 = _bedslopevertices->at(v2index); |
---|
178 | v3 = _bedslopevertices->at(v3index); |
---|
179 | |
---|
180 | side1 = v2 - v1; |
---|
181 | side2 = v3 - v2; |
---|
182 | nrm = side1^side2; |
---|
183 | nrm.normalize(); |
---|
184 | |
---|
185 | _bedslopenormals->push_back( nrm ); |
---|
186 | _bedslopecentroids->push_back( (v1+v2+v3)/3.0 ); |
---|
187 | } |
---|
188 | |
---|
189 | // success! |
---|
190 | _valid = true; |
---|
191 | |
---|
192 | } |
---|
193 | |
---|
194 | |
---|
195 | |
---|
196 | SWWReader::~SWWReader() |
---|
197 | { |
---|
198 | _status.push_back( nc_close(_ncid) ); |
---|
199 | delete[]( _px, _py, _pz, _ptime, _pvolumes, _pstage ); |
---|
200 | } |
---|
201 | |
---|
202 | |
---|
203 | |
---|
204 | osg::ref_ptr<osg::Vec3Array> SWWReader::getStageVertexArray(unsigned int index) |
---|
205 | { |
---|
206 | size_t start[2], count[2], iv; |
---|
207 | const ptrdiff_t stride[2] = {1,1}; |
---|
208 | start[0] = index; |
---|
209 | start[1] = 0; |
---|
210 | count[0] = 1; |
---|
211 | count[1] = _npoints; |
---|
212 | |
---|
213 | // stage heights from netcdf file (x and y are same as bedslope) |
---|
214 | int status = nc_get_vars_float (_ncid, _stageid, start, count, stride, _pstage); |
---|
215 | |
---|
216 | // load stage vertex array, scaling and shifting vertices to lie in the unit cube |
---|
217 | _stagevertices = new osg::Vec3Array; |
---|
218 | for (iv=0; iv < _npoints; iv++) |
---|
219 | _stagevertices->push_back( osg::Vec3( |
---|
220 | (_px[iv]-_xoffset)*_xscale - 0.5, |
---|
221 | (_py[iv]-_yoffset)*_yscale - 0.5, |
---|
222 | (_pstage[iv]-_zoffset)*_zscale - 0.5) ); |
---|
223 | |
---|
224 | // stage index and per primitive normal and centroid arrays |
---|
225 | _stageprimitivenormals = new osg::Vec3Array; |
---|
226 | _stagecentroids = new osg::Vec3Array; |
---|
227 | _stageindices = new osg::UIntArray; |
---|
228 | osg::Vec3 v1b, v2b, v3b; |
---|
229 | osg::Vec3 v1s, v2s, v3s; |
---|
230 | osg::Vec3 side1, side2, nrm; |
---|
231 | unsigned int v1index, v2index, v3index; |
---|
232 | osg::ref_ptr<osg::UIntArray> culledlist = new osg::UIntArray; |
---|
233 | std::vector<triangle_list> stageconnectivity = std::vector<triangle_list>(_npoints); |
---|
234 | |
---|
235 | // over all stage triangles |
---|
236 | unsigned int nonculled = 0; |
---|
237 | for (iv=0; iv < _nvolumes; iv++){ |
---|
238 | |
---|
239 | v1index = _pvolumes[3*iv+0]; |
---|
240 | v2index = _pvolumes[3*iv+1]; |
---|
241 | v3index = _pvolumes[3*iv+2]; |
---|
242 | |
---|
243 | v1s = _stagevertices->at(v1index); |
---|
244 | v2s = _stagevertices->at(v2index); |
---|
245 | v3s = _stagevertices->at(v3index); |
---|
246 | v1b = _bedslopevertices->at(v1index); |
---|
247 | v2b = _bedslopevertices->at(v2index); |
---|
248 | v3b = _bedslopevertices->at(v3index); |
---|
249 | |
---|
250 | // shallow water depth culling test |
---|
251 | if( (_cullsetting == CULLALL || _cullsetting == CULLNEARZERO) && |
---|
252 | ((v1s.z()+v2s.z()+v3s.z())/3.0 - _bedslopecentroids->at(iv).z()) < _cullnearzero ) |
---|
253 | { |
---|
254 | culledlist->push_back( iv ); |
---|
255 | continue; |
---|
256 | } |
---|
257 | |
---|
258 | // steep bedslope culling test |
---|
259 | if( _cullsetting == CULLALL || _cullsetting == CULLSTEEPANGLE ) |
---|
260 | { |
---|
261 | if( _bedslopenormals->at(iv) * osg::Vec3f(0,0,1) < _cullsteepangle ) |
---|
262 | { |
---|
263 | culledlist->push_back( iv ); |
---|
264 | continue; |
---|
265 | } |
---|
266 | } |
---|
267 | |
---|
268 | // pass through, current triangle is visible (not culled) |
---|
269 | _stageindices->push_back( v1index ); |
---|
270 | _stageindices->push_back( v2index ); |
---|
271 | _stageindices->push_back( v3index ); |
---|
272 | |
---|
273 | // current triangle primitive normal |
---|
274 | side1 = v2s - v1s; |
---|
275 | side2 = v3s - v2s; |
---|
276 | nrm = side1^side2; |
---|
277 | nrm.normalize(); |
---|
278 | |
---|
279 | // stage triangle connectivity, a list (indexed by vertex number) |
---|
280 | // of lists (indices of non-culled triangles sharing this vertex) |
---|
281 | stageconnectivity.at(v1index).push_back(nonculled); |
---|
282 | stageconnectivity.at(v2index).push_back(nonculled); |
---|
283 | stageconnectivity.at(v3index).push_back(nonculled); |
---|
284 | |
---|
285 | // store centroid and primitive normal (both needed to visualize vectors) |
---|
286 | _stagecentroids->push_back( (v1s+v2s+v3s)/3.0 ); |
---|
287 | _stageprimitivenormals->push_back( nrm ); |
---|
288 | |
---|
289 | // non-culled triangle count |
---|
290 | nonculled++; |
---|
291 | |
---|
292 | } |
---|
293 | |
---|
294 | // per-vertex normals |
---|
295 | _stagevertexnormals = new osg::Vec3Array; |
---|
296 | int num_shared_triangles, triangle_index; |
---|
297 | for (iv=0; iv < _npoints; iv++) |
---|
298 | { |
---|
299 | nrm.set(0,0,0); |
---|
300 | num_shared_triangles = stageconnectivity.at(iv).size(); |
---|
301 | for (int i=0; i<num_shared_triangles; i++ ) |
---|
302 | { |
---|
303 | triangle_index = stageconnectivity.at(iv).at(i); |
---|
304 | |
---|
305 | // summing of contributing primitive normals |
---|
306 | nrm += _stageprimitivenormals->at(triangle_index); |
---|
307 | } |
---|
308 | |
---|
309 | // FIXME: Length of this list has to match vertex list length. |
---|
310 | // We have vertices and normals that aren't being referenced. |
---|
311 | if( 1 || num_shared_triangles ) // avoid vertices belonging only to culled triangles |
---|
312 | { |
---|
313 | nrm = nrm / num_shared_triangles; // average |
---|
314 | nrm.normalize(); |
---|
315 | _stagevertexnormals->push_back(nrm); |
---|
316 | } |
---|
317 | } |
---|
318 | |
---|
319 | return _stagevertices; |
---|
320 | } |
---|
321 | |
---|
322 | |
---|
323 | |
---|
324 | void SWWReader::toggleCullSetting() |
---|
325 | { |
---|
326 | |
---|
327 | switch( _cullsetting ) |
---|
328 | { |
---|
329 | |
---|
330 | case CULLALL: |
---|
331 | _cullsetting = CULLNEARZERO; |
---|
332 | osg::notify(osg::INFO) << "[SWWReader::toggleCullSetting] CULLNEARZERO" << std::endl; |
---|
333 | break; |
---|
334 | |
---|
335 | case CULLNEARZERO: |
---|
336 | _cullsetting = CULLSTEEPANGLE; |
---|
337 | osg::notify(osg::INFO) << "[SWWReader::toggleCullSetting] CULLSTEEPANGLE" << std::endl; |
---|
338 | break; |
---|
339 | |
---|
340 | case CULLSTEEPANGLE: |
---|
341 | _cullsetting = CULLNONE; |
---|
342 | osg::notify(osg::INFO) << "[SWWReader::toggleCullSetting] CULLNONE" << std::endl; |
---|
343 | break; |
---|
344 | |
---|
345 | case CULLNONE: |
---|
346 | _cullsetting = CULLALL; |
---|
347 | osg::notify(osg::INFO) << "[SWWReader::toggleCullSetting] CULLALL" << std::endl; |
---|
348 | break; |
---|
349 | |
---|
350 | } |
---|
351 | } |
---|
352 | |
---|
353 | |
---|
354 | |
---|
355 | bool SWWReader::_statusHasError() |
---|
356 | { |
---|
357 | bool haserror = false; // assume success, trap failure |
---|
358 | |
---|
359 | std::vector<int>::iterator iter; |
---|
360 | for (iter=_status.begin(); iter != _status.end(); iter++) |
---|
361 | { |
---|
362 | if (*iter != NC_NOERR) |
---|
363 | { |
---|
364 | osg::notify(osg::WARN) << nc_strerror(*iter) << std::endl; |
---|
365 | haserror = true; |
---|
366 | nc_close(_ncid); |
---|
367 | break; |
---|
368 | } |
---|
369 | } |
---|
370 | |
---|
371 | // on return start gathering result values afresh |
---|
372 | _status.clear(); |
---|
373 | |
---|
374 | return haserror; |
---|
375 | } |
---|