source: Swollen/swwreader/SWWReader.cpp @ 6

Last change on this file since 6 was 6, checked in by darran, 20 years ago

new import

File size: 10.0 KB
Line 
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// only constructor, requires netcdf file
38SWWReader::SWWReader(const std::string& filename)
39{
40    // assume failure until end of constructor
41    _valid = false;
42
43    // netcdf filename
44    _filename = new std::string(filename);
45
46    // netcdf open
47    _status.push_back( nc_open(_filename->c_str(), NC_NOWRITE, &_ncid) );
48    if (this->_statusHasError()) return;
49
50    // dimension ids
51    _status.push_back( nc_inq_dimid(_ncid, "number_of_volumes", &_nvolumesid) );
52    _status.push_back( nc_inq_dimid(_ncid, "number_of_vertices", &_nverticesid) );
53    _status.push_back( nc_inq_dimid(_ncid, "number_of_points", &_npointsid) );
54    _status.push_back( nc_inq_dimid(_ncid, "number_of_timesteps", &_ntimestepsid) );
55    if (this->_statusHasError()) return;
56
57    // dimension values
58    _status.push_back( nc_inq_dimlen(_ncid, _nvolumesid, &_nvolumes) );
59    _status.push_back( nc_inq_dimlen(_ncid, _nverticesid, &_nvertices) );
60    _status.push_back( nc_inq_dimlen(_ncid, _npointsid, &_npoints) );
61    _status.push_back( nc_inq_dimlen(_ncid, _ntimestepsid, &_ntimesteps) );
62    if (this->_statusHasError()) return;
63
64    // debug info
65    osg::notify(osg::INFO) << "[SWWReader]" <<  std::endl;
66    osg::notify(osg::INFO) << "    nvolumes: " << _nvolumes <<  std::endl;
67    osg::notify(osg::INFO) << "    nvertices: " << _nvertices << std::endl;
68    osg::notify(osg::INFO) << "    npoints: " << _npoints << std::endl;
69    osg::notify(osg::INFO) << "    ntimesteps: " << _ntimesteps << std::endl;
70
71
72    // variable ids
73    _status.push_back( nc_inq_varid (_ncid, "x", &_xid) );
74    _status.push_back( nc_inq_varid (_ncid, "y", &_yid) );
75    _status.push_back( nc_inq_varid (_ncid, "z", &_zid) );
76    _status.push_back( nc_inq_varid (_ncid, "volumes", &_volumesid) );
77    _status.push_back( nc_inq_varid (_ncid, "time", &_timeid) );
78    _status.push_back( nc_inq_varid (_ncid, "stage", &_stageid) );
79    if (this->_statusHasError()) return;
80
81    // allocation of variable arrays, destructor responsible for cleanup
82    _px = new float[_npoints];
83    _py = new float[_npoints];
84    _pz = new float[_npoints];
85    _ptime = new float[_ntimesteps];
86    _pvolumes = new unsigned int[_nvertices * _nvolumes];
87    _pstage = new float[_npoints];
88
89    // loading variables from netcdf file
90    _status.push_back( nc_get_var_float (_ncid, _xid, _px) );  // x vertices
91    _status.push_back( nc_get_var_float (_ncid, _yid, _py) );  // y vertices
92    _status.push_back( nc_get_var_float (_ncid, _zid, _pz) );  // bedslope heights
93    _status.push_back( nc_get_var_float (_ncid, _timeid, _ptime) );  // time array
94    _status.push_back( nc_get_var_int (_ncid, _volumesid, (int *) _pvolumes) );  // triangle indices
95    if (this->_statusHasError()) return;
96
97    // loop index
98    size_t iv;
99
100    // vertex indices
101    unsigned int v1index, v2index, v3index;
102
103    // culling defaults
104    _cullnearzero = true;
105    _cullnearzerovalue = 0.001;
106    _cullsteepangle = true;
107    _cullsteepanglevalue = 0.9;
108
109    // compute triangle connectivity, a list (indexed by vertex number)
110    // of lists (indices of triangles sharing this vertex)
111    _connectivity = std::vector<triangle_list>(_npoints);
112    for (iv=0; iv < _nvolumes; iv++)
113    {
114        v1index = _pvolumes[3*iv+0];
115        v2index = _pvolumes[3*iv+1];
116        v3index = _pvolumes[3*iv+2];
117        _connectivity.at(v1index).push_back(iv);
118        _connectivity.at(v2index).push_back(iv);
119        _connectivity.at(v3index).push_back(iv);
120    }
121
122    // load bedslope vertex array
123    _bedslopevertices = new osg::Vec3Array;
124    for (iv=0; iv < _npoints; iv++)
125        _bedslopevertices->push_back( osg::Vec3(_px[iv], _py[iv], _pz[iv]) );
126
127    // load bedslope index array, pvolumes array indexes into x, y and z
128    _bedslopeindices = new osg::UIntArray;
129    for (iv=0; iv < _nvolumes*_nvertices; iv++)
130        _bedslopeindices->push_back( _pvolumes[iv] );
131
132    // calculate bedslope primitive normal and centroid arrays
133    _bedslopenormals = new osg::Vec3Array;
134    _bedslopecentroids = new osg::Vec3Array;
135    osg::Vec3 v1, v2, v3, side1, side2, nrm;
136    for (iv=0; iv < _nvolumes; iv++)
137    {
138        v1index = _pvolumes[3*iv+0];
139        v2index = _pvolumes[3*iv+1];
140        v3index = _pvolumes[3*iv+2];
141
142        v1 = _bedslopevertices->at(v1index);
143        v2 = _bedslopevertices->at(v2index);
144        v3 = _bedslopevertices->at(v3index);
145
146        side1 = v2 - v1;
147        side2 = v3 - v2;
148        nrm = side1^side2;
149        nrm.normalize();
150
151        _bedslopenormals->push_back( nrm );
152        _bedslopecentroids->push_back( (v1+v2+v3)/3.0 );
153    }
154
155    // success!
156    _valid = true;
157
158}
159
160
161
162SWWReader::~SWWReader()
163{
164    _status.push_back( nc_close(_ncid) );
165    delete[]( _px, _py, _pz, _ptime, _pvolumes, _pstage );
166}
167
168
169
170osg::ref_ptr<osg::Vec3Array> SWWReader::getStageVertexArray(unsigned int index)
171{
172    size_t start[2], count[2], iv;
173    const ptrdiff_t stride[2] = {1,1};
174    start[0] = index;
175    start[1] = 0;
176    count[0] = 1;
177    count[1] = _npoints;
178
179    // stage heights from netcdf file (x and y are same as bedslope)
180    int status = nc_get_vars_float (_ncid, _stageid, start, count, stride, _pstage);
181
182    // load stage vertex array
183    _stagevertices = new osg::Vec3Array;
184    for (iv=0; iv < _npoints; iv++)
185        _stagevertices->push_back( osg::Vec3(_px[iv], _py[iv], _pstage[iv]) );
186
187    // stage index and per primitive normal and centroid arrays
188    _stageprimitivenormals = new osg::Vec3Array;
189    _stagecentroids = new osg::Vec3Array;
190    _stageindices = new osg::UIntArray;
191    osg::Vec3 v1b, v2b, v3b;
192    osg::Vec3 v1s, v2s, v3s;
193    osg::Vec3 side1, side2, nrm;
194    unsigned int v1index, v2index, v3index;
195    osg::ref_ptr<osg::UIntArray> culledlist = new osg::UIntArray;
196    std::vector<triangle_list> stageconnectivity = std::vector<triangle_list>(_npoints);
197
198    // over all stage triangles
199    unsigned int nonculled = 0;
200    for (iv=0; iv < _nvolumes; iv++){
201
202        v1index = _pvolumes[3*iv+0];
203        v2index = _pvolumes[3*iv+1];
204        v3index = _pvolumes[3*iv+2];
205
206        v1s = _stagevertices->at(v1index);
207        v2s = _stagevertices->at(v2index);
208        v3s = _stagevertices->at(v3index);
209        v1b = _bedslopevertices->at(v1index);
210        v2b = _bedslopevertices->at(v2index);
211        v3b = _bedslopevertices->at(v3index);
212
213        // shallow water depth culling test
214        if( _cullnearzero && ((v1s.z()+v2s.z()+v3s.z())/3.0 - _bedslopecentroids->at(iv).z()) < _cullnearzerovalue )
215        {
216            culledlist->push_back( iv );
217            continue;
218        }
219
220        // steep bedslope culling test
221        if( _cullsteepangle )
222        {
223            if( _bedslopenormals->at(iv) * osg::Vec3f(0,0,1) < _cullsteepanglevalue )
224            {
225                culledlist->push_back( iv );
226                continue;
227            }
228        }
229
230        // pass through, current triangle is visible (not culled)
231        _stageindices->push_back( v1index );
232        _stageindices->push_back( v2index );
233        _stageindices->push_back( v3index );
234
235        // current triangle primitive normal
236        side1 = v2s - v1s;
237        side2 = v3s - v2s;
238        nrm = side1^side2;
239        nrm.normalize();
240
241        // stage triangle connectivity, a list (indexed by vertex number)
242        // of lists (indices of non-culled triangles sharing this vertex)
243        stageconnectivity.at(v1index).push_back(nonculled);
244        stageconnectivity.at(v2index).push_back(nonculled);
245        stageconnectivity.at(v3index).push_back(nonculled);
246
247        // store centroid and primitive normal (both needed to visualize vectors)
248        _stagecentroids->push_back( (v1s+v2s+v3s)/3.0 );
249        _stageprimitivenormals->push_back( nrm );
250
251        // non-culled triangle count
252        nonculled++;
253
254    }
255
256    // per-vertex normals
257    _stagevertexnormals = new osg::Vec3Array;
258    int num_shared_triangles, triangle_index;
259    for (iv=0; iv < _npoints; iv++)
260    {
261        nrm.set(0,0,0);
262        num_shared_triangles = stageconnectivity.at(iv).size();
263        for (int i=0; i<num_shared_triangles; i++ )
264        {
265            triangle_index = stageconnectivity.at(iv).at(i);
266
267            // summing of contributing primitive normals
268            nrm += _stageprimitivenormals->at(triangle_index);
269        }
270
271        // FIXME: Length of this list has to match vertex list length.
272        // We have vertices and normals that aren't being referenced.
273        if( 1 || num_shared_triangles )  // avoid vertices belonging only to culled triangles
274        {
275            nrm = nrm / num_shared_triangles;  // average
276            nrm.normalize();
277            _stagevertexnormals->push_back(nrm);
278        }
279    }
280
281    return _stagevertices;
282}
283
284
285bool SWWReader::_statusHasError()
286{
287    bool haserror = false;  // assume success, trap failure
288
289    std::vector<int>::iterator iter;
290    for (iter=_status.begin(); iter != _status.end(); iter++)
291    {
292        if (*iter != NC_NOERR)
293        {
294            osg::notify(osg::WARN) << nc_strerror(*iter) <<  std::endl;
295            haserror = true;
296            nc_close(_ncid);
297            break;
298        }
299    }
300
301    // on return start gathering result values afresh
302    _status.clear();
303
304    return haserror;
305}
Note: See TracBrowser for help on using the repository browser.