source: Swollen/swwreader/swwreader.cpp @ 92

Last change on this file since 92 was 92, checked in by darran, 19 years ago
  • bedslope texture only when asked
  • code cleanup (comment-out PnP gone)
  • bedslope material/stateset only used for texture
  • 'l' keyboard toggle removed
  • code for looking at global attributes of .sww file
File size: 12.9 KB
Line 
1/*
2    SWWReader
3
4    Reader of SWW files from within OpenSceneGraph.
5
6    copyright (C) 2004 Geoscience Australia
7*/
8
9#include <SWWReader.h>
10#include <string>
11#include <fstream>
12#include <netcdf.h>
13#include <osg/Notify>
14
15#include <osgDB/Registry>
16#include <osgDB/ReadFile>
17
18#include <stdlib.h>
19#include <stdio.h>
20#include <gdal_priv.h>
21
22
23// only constructor, requires netcdf file
24SWWReader::SWWReader(const std::string& filename)
25{
26    // assume failure until end of constructor
27    _valid = false;
28
29    // netcdf filename
30    _filename = new std::string(filename);
31
32    // netcdf open
33    _status.push_back( nc_open(_filename->c_str(), NC_NOWRITE, &_ncid) );
34    if (this->_statusHasError()) return;
35
36    // dimension ids
37    _status.push_back( nc_inq_dimid(_ncid, "number_of_volumes", &_nvolumesid) );
38    _status.push_back( nc_inq_dimid(_ncid, "number_of_vertices", &_nverticesid) );
39    _status.push_back( nc_inq_dimid(_ncid, "number_of_points", &_npointsid) );
40    _status.push_back( nc_inq_dimid(_ncid, "number_of_timesteps", &_ntimestepsid) );
41    if (this->_statusHasError()) return;
42
43    // dimension values
44    _status.push_back( nc_inq_dimlen(_ncid, _nvolumesid, &_nvolumes) );
45    _status.push_back( nc_inq_dimlen(_ncid, _nverticesid, &_nvertices) );
46    _status.push_back( nc_inq_dimlen(_ncid, _npointsid, &_npoints) );
47    _status.push_back( nc_inq_dimlen(_ncid, _ntimestepsid, &_ntimesteps) );
48    if (this->_statusHasError()) return;
49
50    // variable ids
51    _status.push_back( nc_inq_varid (_ncid, "x", &_xid) );
52    _status.push_back( nc_inq_varid (_ncid, "y", &_yid) );
53    _status.push_back( nc_inq_varid (_ncid, "z", &_zid) );
54    _status.push_back( nc_inq_varid (_ncid, "volumes", &_volumesid) );
55    _status.push_back( nc_inq_varid (_ncid, "time", &_timeid) );
56    _status.push_back( nc_inq_varid (_ncid, "stage", &_stageid) );
57    if (this->_statusHasError()) return;
58
59    // allocation of variable arrays, destructor responsible for cleanup
60    _px = new float[_npoints];
61    _py = new float[_npoints];
62    _pz = new float[_npoints];
63    _ptime = new float[_ntimesteps];
64    _pvolumes = new unsigned int[_nvertices * _nvolumes];
65    _pstage = new float[_npoints];
66
67    // loading variables from netcdf file
68    _status.push_back( nc_get_var_float (_ncid, _xid, _px) );  // x vertices
69    _status.push_back( nc_get_var_float (_ncid, _yid, _py) );  // y vertices
70    _status.push_back( nc_get_var_float (_ncid, _zid, _pz) );  // bedslope heights
71    _status.push_back( nc_get_var_float (_ncid, _timeid, _ptime) );  // time array
72    _status.push_back( nc_get_var_int (_ncid, _volumesid, (int *) _pvolumes) );  // triangle indices
73    if (this->_statusHasError()) return;
74
75
76        // sww file can optionally contain bedslope texture image filename
77        size_t attlen;
78        if( nc_inq_attlen(_ncid, NC_GLOBAL, "institution", &attlen) != NC_ENOTATT )
79        {
80                std::string filename;
81                int nc_get_att_text(_ncid, NC_GLOBAL, "institution", char filename);
82                setBedslopeTexture( std::string filename );
83        }
84
85    // loop index
86    size_t iv;
87
88    // vertex indices
89    unsigned int v1index, v2index, v3index;
90
91    // bedslope extents and resultant scale factors
92    float xmin, xmax, xrange;
93    float ymin, ymax, yrange;
94    float zmin, zmax, zrange;
95    float aspect_ratio;
96    xmin = _px[0]; xmax = _px[0];
97    ymin = _py[0]; ymax = _py[0];
98    zmin = _pz[0]; zmax = _pz[0];
99    for( iv=1; iv < _npoints; iv++ )
100    {
101        if( _px[iv] < xmin ) xmin = _px[iv];
102        if( _px[iv] > xmax ) xmax = _px[iv];
103        if( _py[iv] < ymin ) ymin = _py[iv];
104        if( _py[iv] > ymax ) ymax = _py[iv];
105        if( _pz[iv] < zmin ) zmin = _pz[iv];
106        if( _pz[iv] > zmax ) zmax = _pz[iv];
107    }
108
109    xrange = xmax - xmin;
110    yrange = ymax - ymin;
111    zrange = zmax - zmin;
112    aspect_ratio = xrange/yrange;
113    _xscale = 1.0 / xrange;
114    _yscale = 1.0 / yrange;
115   
116    // handle case of a flat bed that doesn't necessarily pass through z=0
117    _zscale = (zrange == 0.0) ? 1.0 : 1.0/zrange;
118
119    _xoffset = xmin;
120    _yoffset = ymin;
121    _zoffset = zmin;
122    _xcenter = 0.5;
123    _ycenter = 0.5;
124    _zcenter = 0.0;
125
126    if( aspect_ratio > 1.0 )
127    {
128        _scale = 1.0 / xrange;
129        _ycenter /= aspect_ratio;
130    }
131    else
132    {
133        _scale = 1.0 / yrange;
134        _xcenter /= aspect_ratio;
135    }
136
137    // alpha-scaling defaults
138    _alphamin = 0.8;       
139    _alphamax = 1.0;
140    _heightmin = 0.0;
141    _heightmax = 1.0;   
142
143    osg::notify(osg::INFO) << "[SWWReader] xmin: " << xmin <<  std::endl;
144    osg::notify(osg::INFO) << "[SWWReader] xmax: " << xmax <<  std::endl;
145    osg::notify(osg::INFO) << "[SWWReader] xscale: " << _xscale <<  std::endl;
146    osg::notify(osg::INFO) << "[SWWReader] ymin: " << ymin <<  std::endl;
147    osg::notify(osg::INFO) << "[SWWReader] ymax: " << ymax <<  std::endl;
148    osg::notify(osg::INFO) << "[SWWReader] yscale: " << _yscale <<  std::endl;
149    osg::notify(osg::INFO) << "[SWWReader] zmin: " << zmin <<  std::endl;
150    osg::notify(osg::INFO) << "[SWWReader] zmax: " << zmax <<  std::endl;
151    osg::notify(osg::INFO) << "[SWWReader] zscale: " << _zscale <<  std::endl;
152
153
154    // compute triangle connectivity, a list (indexed by vertex number)
155    // of lists (indices of triangles sharing this vertex)
156    _connectivity = std::vector<triangle_list>(_npoints);
157    for (iv=0; iv < _nvolumes; iv++)
158    {
159        v1index = _pvolumes[3*iv+0];
160        v2index = _pvolumes[3*iv+1];
161        v3index = _pvolumes[3*iv+2];
162        _connectivity.at(v1index).push_back(iv);
163        _connectivity.at(v2index).push_back(iv);
164        _connectivity.at(v3index).push_back(iv);
165    }
166
167    // bedslope vertex array, shifting and scaling vertices to unit cube
168    // centred about the origin
169    _bedslopevertices = new osg::Vec3Array;
170    for (iv=0; iv < _npoints; iv++)
171        _bedslopevertices->push_back( osg::Vec3( (_px[iv]-_xoffset)*_scale - _xcenter, 
172                                                 (_py[iv]-_yoffset)*_scale - _ycenter, 
173                                                 (_pz[iv]-_zoffset)*_scale - _zcenter) );
174
175    // bedslope index array, pvolumes array indexes into x, y and z
176    _bedslopeindices = new osg::DrawElementsUShort(osg::PrimitiveSet::TRIANGLES, _nvolumes*_nvertices);
177    for (iv=0; iv < _nvolumes*_nvertices; iv++)
178        (*_bedslopeindices)[iv] = _pvolumes[iv];
179
180
181    // calculate bedslope primitive normal and centroid arrays
182    _bedslopenormals = new osg::Vec3Array;
183    _bedslopecentroids = new osg::Vec3Array;
184    osg::Vec3 v1, v2, v3, side1, side2, nrm;
185    for (iv=0; iv < _nvolumes; iv++)
186    {
187        v1index = _pvolumes[3*iv+0];
188        v2index = _pvolumes[3*iv+1];
189        v3index = _pvolumes[3*iv+2];
190
191        v1 = _bedslopevertices->at(v1index);
192        v2 = _bedslopevertices->at(v2index);
193        v3 = _bedslopevertices->at(v3index);
194
195        side1 = v2 - v1;
196        side2 = v3 - v2;
197        nrm = side1^side2;
198        nrm.normalize();
199
200        _bedslopenormals->push_back( nrm );
201        _bedslopecentroids->push_back( (v1+v2+v3)/3.0 );
202    }
203
204    // success!
205    _valid = true;
206
207}
208
209
210
211SWWReader::~SWWReader()
212{
213    _status.push_back( nc_close(_ncid) );
214    delete[]( _px, _py, _pz, _ptime, _pvolumes, _pstage );
215}
216
217
218
219void SWWReader::setBedslopeTexture( std::string filename )
220{
221    osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] filename: " << filename <<  std::endl;
222   _bedslopetexture = new std::string(filename);
223   _bedslopegeodata.hasData = false;
224
225   // GDAL Geospatial Image Library
226   GDALDataset *pGDAL;
227   GDALAllRegister();
228   pGDAL = (GDALDataset *) GDALOpen( filename.c_str(), GA_ReadOnly );
229   if( pGDAL == NULL )
230      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] GDAL unable to open file: " << filename <<  std::endl;
231   else
232   {
233      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] GDAL Driver: "
234                             << pGDAL->GetDriver()->GetDescription() << "/" 
235                             << pGDAL->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) 
236                             << std::endl;
237
238      _bedslopegeodata.xresolution = pGDAL->GetRasterXSize();
239      _bedslopegeodata.yresolution = pGDAL->GetRasterYSize();
240
241      // check if image contains geodata
242      double geodata[6];
243      if( pGDAL->GetGeoTransform( geodata ) == CE_None )
244      {
245          _bedslopegeodata.xorigin = geodata[0];
246          _bedslopegeodata.yorigin = geodata[3];
247          _bedslopegeodata.rotation = geodata[2];
248          _bedslopegeodata.xpixel = geodata[1];
249          _bedslopegeodata.ypixel = geodata[5];
250          _bedslopegeodata.hasData = true;
251
252          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xorigin: " << _bedslopegeodata.xorigin <<  std::endl;
253          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] yorigin: " << _bedslopegeodata.yorigin <<  std::endl;
254          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xpixel: " << _bedslopegeodata.xpixel <<  std::endl;
255          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] ypixel: " << _bedslopegeodata.ypixel <<  std::endl;
256          osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] rotation: " << _bedslopegeodata.rotation <<  std::endl;
257      }
258   }
259
260}
261
262
263
264osg::Image* SWWReader::getBedslopeTexture()
265{
266   return osgDB::readImageFile( _bedslopetexture->c_str() );
267}
268
269
270
271osg::ref_ptr<osg::Vec2Array> SWWReader::getBedslopeTextureCoords()
272{
273
274   // bedslope texture coords based on (x,y) locations scaled by extents into range [0,1]
275   _bedslopetexcoords = new osg::Vec2Array;
276   for( unsigned int iv=0; iv < _npoints; iv++ )
277      _bedslopetexcoords->push_back( osg::Vec2( (_px[iv]-_xoffset)*_xscale, (_py[iv]-_yoffset)*_yscale ) );
278
279   return _bedslopetexcoords;
280}
281
282
283
284osg::ref_ptr<osg::Vec3Array> SWWReader::getStageVertexArray(unsigned int index)
285{
286    size_t start[2], count[2], iv;
287    const ptrdiff_t stride[2] = {1,1};
288    start[0] = index;
289    start[1] = 0;
290    count[0] = 1;
291    count[1] = _npoints;
292
293    // stage heights from netcdf file (x and y are same as bedslope)
294    int status = nc_get_vars_float (_ncid, _stageid, start, count, stride, _pstage);
295
296    // load stage vertex array, scaling and shifting vertices to lie in the unit cube
297    _stagevertices = new osg::Vec3Array;
298    for (iv=0; iv < _npoints; iv++)
299        _stagevertices->push_back( osg::Vec3( (_px[iv]-_xoffset)*_scale - _xcenter, 
300                                              (_py[iv]-_yoffset)*_scale - _ycenter, 
301                                              (_pstage[iv]-_zoffset)*_scale - _zcenter) );
302
303    // stage index, per primitive normal and centroid arrays
304    _stageprimitivenormals = new osg::Vec3Array;
305    _stagecentroids = new osg::Vec3Array;
306    osg::Vec3 v1b, v2b, v3b;
307    osg::Vec3 v1s, v2s, v3s;
308    osg::Vec3 side1, side2, nrm;
309    unsigned int v1index, v2index, v3index;
310
311    // over all stage triangles
312    for (iv=0; iv < _nvolumes; iv++)
313    {
314        v1index = _pvolumes[3*iv+0];
315        v2index = _pvolumes[3*iv+1];
316        v3index = _pvolumes[3*iv+2];
317
318        v1s = _stagevertices->at(v1index);
319        v2s = _stagevertices->at(v2index);
320        v3s = _stagevertices->at(v3index);
321        v1b = _bedslopevertices->at(v1index);
322        v2b = _bedslopevertices->at(v2index);
323        v3b = _bedslopevertices->at(v3index);
324
325        // current triangle primitive normal
326        side1 = v2s - v1s;
327        side2 = v3s - v2s;
328        nrm = side1^side2;
329        nrm.normalize();
330
331        // store centroid and primitive normal
332        _stagecentroids->push_back( (v1s+v2s+v3s)/3.0 );
333        _stageprimitivenormals->push_back( nrm );
334    }
335
336
337    // stage height above bedslope mapped as alpha value
338    //      alpha = min( a(h-hmin) + alphamin, alphamax),  h >= hmin
339    //      alpha = 0,                                     h < hmin
340    // where a = (alphamax-alphamin)/(hmax-hmin)
341    float alpha, height, alphascale;
342    alphascale = (_alphamax - _alphamin) / (_heightmax - _heightmin);
343    _stagecolors = new osg::Vec4Array;
344    for (iv=0; iv < _npoints; iv++)
345    {
346        height = _stagevertices->at(iv).z() - _bedslopevertices->at(iv).z();
347       
348        if (height < _heightmin)
349        {
350            alpha = 0.0;
351        } 
352        else 
353        {
354            alpha = alphascale * (height - _heightmin) + _alphamin;
355            if( alpha > _alphamax ) 
356                alpha = _alphamax;       
357        }
358        _stagecolors->push_back( osg::Vec4( 1.0, 1.0, 1.0, alpha ) );
359    }
360
361
362    // per-vertex normals calculated as average of primitive normals
363    // from contributing triangles
364    _stagevertexnormals = new osg::Vec3Array;
365    int num_shared_triangles, triangle_index;
366    for (iv=0; iv < _npoints; iv++)
367    {
368        nrm.set(0,0,0);
369        num_shared_triangles = _connectivity.at(iv).size();
370        for (int i=0; i<num_shared_triangles; i++ )
371        {
372            triangle_index = _connectivity.at(iv).at(i);
373            nrm += _stageprimitivenormals->at(triangle_index);
374        }
375
376        nrm = nrm / num_shared_triangles;  // average
377        nrm.normalize();
378        _stagevertexnormals->push_back(nrm);
379    }
380
381    return _stagevertices;
382}
383
384
385
386
387bool SWWReader::_statusHasError()
388{
389    bool haserror = false;  // assume success, trap failure
390
391    std::vector<int>::iterator iter;
392    for (iter=_status.begin(); iter != _status.end(); iter++)
393    {
394        if (*iter != NC_NOERR)
395        {
396            osg::notify(osg::WARN) << nc_strerror(*iter) <<  std::endl;
397            haserror = true;
398            nc_close(_ncid);
399            break;
400        }
401    }
402
403    // on return start gathering result values afresh
404    _status.clear();
405
406    return haserror;
407}
Note: See TracBrowser for help on using the repository browser.