source: Swollen/swwreader/swwreader.cpp @ 111

Last change on this file since 111 was 111, checked in by darran, 18 years ago
  • refactoring to facilitate storing state
  • segfaulting at the moment ...
File size: 18.1 KB
Line 
1/*
2  SWWReader
3
4  Reader of SWW files from within OpenSceneGraph.
5
6  copyright (C) 2004-2005 Geoscience Australia
7*/
8
9#include <SWWReader.h>
10#include <cstdlib>
11#include <string>
12#include <fstream>
13#include <netcdf.h>
14#include <osg/Notify>
15
16#include <osgDB/Registry>
17#include <osgDB/ReadFile>
18#include <osgDB/FileNameUtils>
19
20#include <stdlib.h>
21#include <stdio.h>
22#include <gdal_priv.h>
23
24
25// compile time defaults
26#define DEFAULT_CULLANGLE 85.0
27#define DEFAULT_ALPHAMIN 0.8
28#define DEFAULT_ALPHAMAX 1.0
29#define DEFAULT_HEIGHTMIN 0.0
30#define DEFAULT_HEIGHTMAX 1.0
31#define DEFAULT_BEDSLOPEOFFSET 0.0
32#define DEFAULT_CULLONSTART false
33
34
35// only constructor, requires netcdf file
36SWWReader::SWWReader(const std::string& filename)
37{
38   // assume failure until end of constructor
39   _valid = false;
40
41   // netcdf filename
42   _state.swwfilename = new std::string(filename);
43
44   // netcdf open
45   _status.push_back( nc_open(filename.c_str(), NC_NOWRITE, &_ncid) );
46   if (this->_statusHasError()) return;
47
48   // dimension ids
49   _status.push_back( nc_inq_dimid(_ncid, "number_of_volumes", &_nvolumesid) );
50   _status.push_back( nc_inq_dimid(_ncid, "number_of_vertices", &_nverticesid) );
51   _status.push_back( nc_inq_dimid(_ncid, "number_of_points", &_npointsid) );
52   _status.push_back( nc_inq_dimid(_ncid, "number_of_timesteps", &_ntimestepsid) );
53   if (this->_statusHasError()) return;
54
55   // dimension values
56   _status.push_back( nc_inq_dimlen(_ncid, _nvolumesid, &_nvolumes) );
57   _status.push_back( nc_inq_dimlen(_ncid, _nverticesid, &_nvertices) );
58   _status.push_back( nc_inq_dimlen(_ncid, _npointsid, &_npoints) );
59   _status.push_back( nc_inq_dimlen(_ncid, _ntimestepsid, &_ntimesteps) );
60   if (this->_statusHasError()) return;
61
62   // variable ids
63   _status.push_back( nc_inq_varid (_ncid, "x", &_xid) );
64   _status.push_back( nc_inq_varid (_ncid, "y", &_yid) );
65   _status.push_back( nc_inq_varid (_ncid, "z", &_zid) );
66   _status.push_back( nc_inq_varid (_ncid, "volumes", &_volumesid) );
67   _status.push_back( nc_inq_varid (_ncid, "time", &_timeid) );
68   _status.push_back( nc_inq_varid (_ncid, "stage", &_stageid) );
69   if (this->_statusHasError()) return;
70
71   // allocation of variable arrays, destructor responsible for cleanup
72   _px = new float[_npoints];
73   _py = new float[_npoints];
74   _pz = new float[_npoints];
75   _ptime = new float[_ntimesteps];
76   _pvolumes = new unsigned int[_nvertices * _nvolumes];
77   _pstage = new float[_npoints];
78
79   // loading variables from netcdf file
80   _status.push_back( nc_get_var_float (_ncid, _xid, _px) );  // x vertices
81   _status.push_back( nc_get_var_float (_ncid, _yid, _py) );  // y vertices
82   _status.push_back( nc_get_var_float (_ncid, _zid, _pz) );  // bedslope heights
83   _status.push_back( nc_get_var_float (_ncid, _timeid, _ptime) );  // time array
84   _status.push_back( nc_get_var_int (_ncid, _volumesid, (int *) _pvolumes) );  // triangle indices
85   if (this->_statusHasError()) return;
86
87
88   osg::notify(osg::INFO) << "[SWWReader] number of volumes: " << _nvolumes <<  std::endl;
89   osg::notify(osg::INFO) << "[SWWReader] number of vertices: " << _nvertices <<  std::endl;
90   osg::notify(osg::INFO) << "[SWWReader] number of points: " << _npoints <<  std::endl;
91   osg::notify(osg::INFO) << "[SWWReader] number of timesteps: " << _ntimesteps <<  std::endl;
92
93
94   // sww file can optionally contain bedslope texture image filename
95   size_t attlen; // length of text attribute (if it exists)
96   if( nc_inq_attlen(_ncid, NC_GLOBAL, "texture", &attlen) != NC_ENOTATT )
97   {
98      char* texfilename;
99      if( (texfilename = (char*) malloc(attlen+1)) != NULL){
100         int status;
101         status = nc_get_att_text(_ncid, NC_GLOBAL, "texture", texfilename);
102         if( status == NC_NOERR )
103         {
104            texfilename[attlen] = '\0';  // ensure string is terminated, not a requirement for netcdf attributes
105            osg::notify(osg::INFO) << "[SWWReader] embedded image filename: " << texfilename <<  std::endl;
106
107            // if sww isn't in current directory, need to prepend sww path to the bedslope texture
108            if( osgDB::getFilePath(filename) == "" )
109               setBedslopeTexture( std::string(texfilename) );
110            else
111               setBedslopeTexture( osgDB::getFilePath(filename) + std::string("/") + std::string(texfilename) );
112         }
113         free( texfilename );
114      }
115   }
116
117   // sww file can optionally contain georeference offset
118   int status1 = nc_get_att_float(_ncid, NC_GLOBAL, "xllcorner", &_xllcorner);
119   int status2 = nc_get_att_float(_ncid, NC_GLOBAL, "yllcorner", &_yllcorner);
120   if( status1 == NC_NOERR && status2 == NC_NOERR)
121   {
122      osg::notify(osg::INFO) << "[SWWReader] xllcorner: " << _xllcorner <<  std::endl;
123      osg::notify(osg::INFO) << "[SWWReader] yllcorner: " << _yllcorner <<  std::endl;
124   }
125   else
126   {
127      _xllcorner = 0.0;  // default value
128      _yllcorner = 0.0;
129   }
130
131
132   // alpha-scaling defaults, can be overridden after construction by command line parameters
133   _state.alphamin = DEFAULT_ALPHAMIN;
134   _state.alphamax = DEFAULT_ALPHAMAX;
135   _state.heightmin = DEFAULT_HEIGHTMIN;
136   _state.heightmax = DEFAULT_HEIGHTMAX;
137
138   // steepness culling default, can be overridden after construction by command line parameter
139   _state.cullangle = DEFAULT_CULLANGLE;
140   _state.culling = DEFAULT_CULLONSTART;
141
142   // loop index
143   size_t iv;
144
145   // vertex indices
146   unsigned int v1index, v2index, v3index;
147
148   // bedslope range and resultant scale factors (note, these don't take into
149   // account any xllcorner, yllcorner offsets - used during texture coord assignment)
150   float xmin, xmax, xrange;
151   float ymin, ymax, yrange;
152   float zmin, zmax, zrange;
153   float aspect_ratio;
154   xmin = _px[0]; 
155   xmax = _px[0];
156   ymin = _py[0]; 
157   ymax = _py[0];
158   zmin = _pz[0]; 
159   zmax = _pz[0];
160   for( iv=1; iv < _npoints; iv++ )
161   {
162      if( _px[iv] < xmin ) xmin = _px[iv];
163      if( _px[iv] > xmax ) xmax = _px[iv];
164      if( _py[iv] < ymin ) ymin = _py[iv];
165      if( _py[iv] > ymax ) ymax = _py[iv];
166      if( _pz[iv] < zmin ) zmin = _pz[iv];
167      if( _pz[iv] > zmax ) zmax = _pz[iv];
168   }
169
170   xrange = xmax - xmin;
171   yrange = ymax - ymin;
172   zrange = zmax - zmin;
173   aspect_ratio = xrange/yrange;
174   _xscale = 1.0 / xrange;
175   _yscale = 1.0 / yrange;
176   
177   // handle case of a flat bed that doesn't necessarily pass through z=0
178   _zscale = (zrange == 0.0) ? 1.0 : 1.0/zrange;
179
180   _xoffset = xmin;
181   _yoffset = ymin;
182   _zoffset = zmin;
183   _xcenter = 0.5;
184   _ycenter = 0.5;
185   _zcenter = 0.0;
186
187   if( aspect_ratio > 1.0 )
188   {
189      _scale = 1.0/xrange;
190      _ycenter /= aspect_ratio;
191   }
192   else
193   {
194      _scale = 1.0/yrange;
195      _xcenter *= aspect_ratio;
196   }
197
198   osg::notify(osg::INFO) << "[SWWReader] xmin: " << xmin <<  std::endl;
199   osg::notify(osg::INFO) << "[SWWReader] xmax: " << xmax <<  std::endl;
200   osg::notify(osg::INFO) << "[SWWReader] xrange: " << xrange <<  std::endl;
201   osg::notify(osg::INFO) << "[SWWReader] xscale: " << _xscale <<  std::endl;
202   osg::notify(osg::INFO) << "[SWWReader] xcenter: " << _xcenter <<  std::endl;
203   osg::notify(osg::INFO) <<  std::endl;
204
205   osg::notify(osg::INFO) << "[SWWReader] ymin: " << ymin <<  std::endl;
206   osg::notify(osg::INFO) << "[SWWReader] ymax: " << ymax <<  std::endl;
207   osg::notify(osg::INFO) << "[SWWReader] yrange: " << yrange <<  std::endl;
208   osg::notify(osg::INFO) << "[SWWReader] yscale: " << _yscale <<  std::endl;
209   osg::notify(osg::INFO) << "[SWWReader] ycenter: " << _ycenter <<  std::endl;
210   osg::notify(osg::INFO) <<  std::endl;
211
212   osg::notify(osg::INFO) << "[SWWReader] zmin: " << zmin <<  std::endl;
213   osg::notify(osg::INFO) << "[SWWReader] zmax: " << zmax <<  std::endl;
214   osg::notify(osg::INFO) << "[SWWReader] zrange: " << zrange <<  std::endl;
215   osg::notify(osg::INFO) << "[SWWReader] zscale: " << _zscale <<  std::endl;
216   osg::notify(osg::INFO) << "[SWWReader] zcenter: " << _zcenter <<  std::endl;
217   osg::notify(osg::INFO) <<  std::endl;
218
219   // compute triangle connectivity, a list (indexed by vertex number)
220   // of lists (indices of triangles sharing this vertex)
221   _connectivity = std::vector<triangle_list>(_npoints);
222   for (iv=0; iv < _nvolumes; iv++)
223   {
224      v1index = _pvolumes[3*iv+0];
225      v2index = _pvolumes[3*iv+1];
226      v3index = _pvolumes[3*iv+2];
227      _connectivity.at(v1index).push_back(iv);
228      _connectivity.at(v2index).push_back(iv);
229      _connectivity.at(v3index).push_back(iv);
230   }
231
232   // bedslope vertex array, shifting and scaling vertices to unit cube
233   // centred about the origin
234   _bedslopevertices = new osg::Vec3Array;
235   for (iv=0; iv < _npoints; iv++)
236      _bedslopevertices->push_back( osg::Vec3( (_px[iv]-_xoffset)*_scale - _xcenter, 
237                                               (_py[iv]-_yoffset)*_scale - _ycenter, 
238                                               (_pz[iv]-_zoffset)*_scale - _zcenter - DEFAULT_BEDSLOPEOFFSET) );
239
240   // bedslope index array, pvolumes array indexes into x, y and z
241   _bedslopeindices = new osg::DrawElementsUInt(osg::PrimitiveSet::TRIANGLES, _nvolumes*_nvertices);
242   for (iv=0; iv < _nvolumes*_nvertices; iv++)
243      (*_bedslopeindices)[iv] = _pvolumes[iv];
244
245
246   // calculate bedslope primitive normal and centroid arrays
247   _bedslopenormals = new osg::Vec3Array;
248   _bedslopecentroids = new osg::Vec3Array;
249   osg::Vec3 v1, v2, v3, side1, side2, nrm;
250   for (iv=0; iv < _nvolumes; iv++)
251   {
252      v1index = _pvolumes[3*iv+0];
253      v2index = _pvolumes[3*iv+1];
254      v3index = _pvolumes[3*iv+2];
255
256      v1 = _bedslopevertices->at(v1index);
257      v2 = _bedslopevertices->at(v2index);
258      v3 = _bedslopevertices->at(v3index);
259
260      side1 = v2 - v1;
261      side2 = v3 - v2;
262      nrm = side1^side2;
263      nrm.normalize();
264
265      _bedslopenormals->push_back( nrm );
266      _bedslopecentroids->push_back( (v1+v2+v3)/3.0 );
267   }
268
269   // success!
270   _valid = true;
271
272   std::cout << _state.alphamin << std::endl;
273   std::cout << _state.alphamax << std::endl;
274   std::cout << _state.heightmin << std::endl;
275   std::cout << _state.heightmax << std::endl;
276   std::cout << _state.cullangle << std::endl;
277   std::cout << _state.culling << std::endl;
278
279   std::cout << *_state.swwfilename << std::endl;
280   std::cout << *_state.bedslopetexturefilename << std::endl;
281
282}
283
284
285
286SWWReader::~SWWReader()
287{
288   _status.push_back( nc_close(_ncid) );
289   delete[]( _px, _py, _pz, _ptime, _pvolumes, _pstage );
290}
291
292
293
294void SWWReader::setBedslopeTexture( std::string filename )
295{
296   osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] filename: " << filename <<  std::endl;
297   _state.bedslopetexturefilename = new std::string(filename);
298   _bedslopegeodata.hasData = false;
299
300   // GDAL Geospatial Image Library
301   GDALDataset *pGDAL;
302   GDALAllRegister();
303   pGDAL = (GDALDataset *) GDALOpen( filename.c_str(), GA_ReadOnly );
304   if( pGDAL == NULL )
305      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] GDAL unable to open file: " << filename <<  std::endl;
306   else
307   {
308      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] GDAL Driver: "
309                             << pGDAL->GetDriver()->GetDescription() << "/" 
310                             << pGDAL->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) 
311                             << std::endl;
312
313      _bedslopegeodata.xresolution = pGDAL->GetRasterXSize();
314      _bedslopegeodata.yresolution = pGDAL->GetRasterYSize();
315      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xresolution: " << _bedslopegeodata.xresolution <<  std::endl;
316      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] yresolution: " << _bedslopegeodata.yresolution <<  std::endl;
317
318      // check if image contains geodata
319      double geodata[6];
320      if( pGDAL->GetGeoTransform( geodata ) == CE_None )
321      {
322         _bedslopegeodata.xorigin = geodata[0];
323         _bedslopegeodata.yorigin = geodata[3];
324         _bedslopegeodata.rotation = geodata[2];
325         _bedslopegeodata.xpixel = geodata[1];
326         _bedslopegeodata.ypixel = geodata[5];
327         _bedslopegeodata.hasData = true;
328
329         osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xorigin: " << _bedslopegeodata.xorigin <<  std::endl;
330         osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] yorigin: " << _bedslopegeodata.yorigin <<  std::endl;
331         osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xpixel: " << _bedslopegeodata.xpixel <<  std::endl;
332         osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] ypixel: " << _bedslopegeodata.ypixel <<  std::endl;
333         osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] rotation: " << _bedslopegeodata.rotation <<  std::endl;
334      }
335   }
336}
337
338
339
340osg::Image* SWWReader::getBedslopeTexture()
341{
342   return osgDB::readImageFile( _state.bedslopetexturefilename->c_str() );
343}
344
345
346
347osg::ref_ptr<osg::Vec2Array> SWWReader::getBedslopeTextureCoords()
348{
349   _bedslopetexcoords = new osg::Vec2Array;
350   if( _bedslopegeodata.hasData )  // georeferenced bedslope texture
351   {
352      float xorigin = _bedslopegeodata.xorigin;
353      float yorigin = _bedslopegeodata.yorigin;
354      float xrange = _bedslopegeodata.xpixel * _bedslopegeodata.xresolution;
355      float yrange = _bedslopegeodata.ypixel * _bedslopegeodata.yresolution;
356      for( unsigned int iv=0; iv < _npoints; iv++ )
357         _bedslopetexcoords->push_back( osg::Vec2( (_px[iv] + _xllcorner - xorigin)/xrange, 
358                                                   (_py[iv] + _yllcorner - yorigin)/yrange ) );
359   }
360   else
361   {
362      // decaled using (x,y) locations scaled by extents into range [0,1]
363      for( unsigned int iv=0; iv < _npoints; iv++ )
364         _bedslopetexcoords->push_back( osg::Vec2( (_px[iv]-_xoffset)*_xscale, (_py[iv]-_yoffset)*_yscale ) );
365   }
366
367   return _bedslopetexcoords;
368}
369
370
371
372osg::ref_ptr<osg::Vec3Array> SWWReader::getStageVertexArray(unsigned int index)
373{
374   size_t start[2], count[2], iv;
375   const ptrdiff_t stride[2] = {1,1};
376   start[0] = index;
377   start[1] = 0;
378   count[0] = 1;
379   count[1] = _npoints;
380
381   // empty array for storing list of steep triangles
382   osg::ref_ptr<osg::IntArray> steeptri = new osg::IntArray;
383
384   // stage heights from netcdf file (x and y are same as bedslope)
385   int status = nc_get_vars_float (_ncid, _stageid, start, count, stride, _pstage);
386
387   // load stage vertex array, scaling and shifting vertices to lie in the unit cube
388   _stagevertices = new osg::Vec3Array;
389   for (iv=0; iv < _npoints; iv++)
390      _stagevertices->push_back( osg::Vec3( (_px[iv]-_xoffset)*_scale - _xcenter, 
391                                            (_py[iv]-_yoffset)*_scale - _ycenter, 
392                                            (_pstage[iv]-_zoffset)*_scale - _zcenter) );
393
394   // stage index, per primitive normal and centroid arrays
395   _stageprimitivenormals = new osg::Vec3Array;
396   _stagecentroids = new osg::Vec3Array;
397   osg::Vec3 v1b, v2b, v3b;
398   osg::Vec3 v1s, v2s, v3s;
399   osg::Vec3 side1, side2, nrm;
400   unsigned int v1index, v2index, v3index;
401
402   // cullangle given in degrees, test is against dot product
403   float cullthreshold = cos(osg::DegreesToRadians(_state.cullangle));
404
405   // over all stage triangles
406   for (iv=0; iv < _nvolumes; iv++)
407   {
408      v1index = _pvolumes[3*iv+0];
409      v2index = _pvolumes[3*iv+1];
410      v3index = _pvolumes[3*iv+2];
411
412      v1s = _stagevertices->at(v1index);
413      v2s = _stagevertices->at(v2index);
414      v3s = _stagevertices->at(v3index);
415      v1b = _bedslopevertices->at(v1index);
416      v2b = _bedslopevertices->at(v2index);
417      v3b = _bedslopevertices->at(v3index);
418
419      // current triangle primitive normal
420      side1 = v2s - v1s;
421      side2 = v3s - v2s;
422      nrm = side1^side2;
423      nrm.normalize();
424
425      // store centroid and primitive normal
426      _stagecentroids->push_back( (v1s+v2s+v3s)/3.0 );
427      _stageprimitivenormals->push_back( nrm );
428
429      // identify steep triangles, store index
430      if( fabs(nrm * osg::Vec3f(0,0,1)) < cullthreshold )
431         steeptri->push_back( iv );
432   }
433
434   // stage height above bedslope mapped as alpha value
435   //      alpha = min( a(h-hmin) + alphamin, alphamax),  h >= hmin
436   //      alpha = 0,                                     h < hmin
437   // where a = (alphamax-alphamin)/(hmax-hmin)
438   float alpha, height, alphascale;
439   alphascale = (_state.alphamax - _state.alphamin) / (_state.heightmax - _state.heightmin);
440   _stagecolors = new osg::Vec4Array;
441   for (iv=0; iv < _npoints; iv++)
442   {
443      // water height above corresponding bedslope
444      height = _stagevertices->at(iv).z() - _bedslopevertices->at(iv).z();
445       
446      if (height < _state.heightmin)
447         alpha = 0.0;
448      else 
449      {
450         alpha = alphascale * (height - _state.heightmin) + _state.alphamin;
451         if( alpha > _state.alphamax ) 
452            alpha = _state.alphamax;
453      }
454
455      _stagecolors->push_back( osg::Vec4( 1.0, 1.0, 1.0, alpha ) );
456   }
457
458   // steep triangle vertices should have alpha=0, overwrite such vertex colours
459   if( _state.culling )
460   {
461      for (iv=0; iv < steeptri->size() ; iv++)
462      {
463         int tri = steeptri->at(iv);
464         v1index = _pvolumes[3*tri+0];
465         v2index = _pvolumes[3*tri+1];
466         v3index = _pvolumes[3*tri+2];
467
468         _stagecolors->at(v1index) = osg::Vec4( 1, 1, 1, 0 );
469         _stagecolors->at(v2index) = osg::Vec4( 1, 1, 1, 0 );
470         _stagecolors->at(v3index) = osg::Vec4( 1, 1, 1, 0 );
471      }
472   }
473
474   // per-vertex normals calculated as average of primitive normals
475   // from contributing triangles
476   _stagevertexnormals = new osg::Vec3Array;
477   int num_shared_triangles, triangle_index;
478   for (iv=0; iv < _npoints; iv++)
479   {
480      nrm.set(0,0,0);
481      num_shared_triangles = _connectivity.at(iv).size();
482      for (int i=0; i<num_shared_triangles; i++ )
483      {
484         triangle_index = _connectivity.at(iv).at(i);
485         nrm += _stageprimitivenormals->at(triangle_index);
486      }
487
488      nrm = nrm / num_shared_triangles;  // average
489      nrm.normalize();
490      _stagevertexnormals->push_back(nrm);
491   }
492
493   return _stagevertices;
494}
495
496
497
498
499bool SWWReader::_statusHasError()
500{
501   bool haserror = false;  // assume success, trap failure
502
503   std::vector<int>::iterator iter;
504   for (iter=_status.begin(); iter != _status.end(); iter++)
505   {
506      if (*iter != NC_NOERR)
507      {
508         osg::notify(osg::WARN) << nc_strerror(*iter) <<  std::endl;
509         haserror = true;
510         nc_close(_ncid);
511         break;
512      }
513   }
514
515   // on return start gathering result values afresh
516   _status.clear();
517
518   return haserror;
519}
Note: See TracBrowser for help on using the repository browser.