source: Swollen/swwreader/swwreader.cpp @ 9

/*
    SWWReader

    Reader of SWW files from within OpenSceneGraph.

    An SWW file is the visualization output of the pyVolution
    Shallow Water Wave equation solver.  SWW files contain bedslope
    geometry and a sequence of stages at known timesteps.  The
    internal format is NetCDF:

    netcdf demo.sww {
        dimensions:
            number_of_volumes = <integer> ;
            number_of_vertices = 3 ;
            number_of_points = <integer> ;
            number_of_timesteps = UNLIMITED ;
        variables:
            float x(number_of_points) ;
            float y(number_of_points) ;
            float z(number_of_points) ;
            int volumes(number_of_volumes, number_of_vertices) ;
            float time(number_of_timesteps) ;
            float stage(number_of_timesteps, number_of_points) ;

    copyright (C) 2004 Geoscience Australia
*/

#include <SWWReader.h>
#include <string>
#include <fstream>
#include <netcdf.h>
#include <osg/Notify>

#include <stdlib.h>
#include <stdio.h>

// only constructor, requires netcdf file
SWWReader::SWWReader(const std::string& filename)
{
    // assume failure until end of constructor
    _valid = false;

    // netcdf filename
    _filename = new std::string(filename);

    // netcdf open
    _status.push_back( nc_open(_filename->c_str(), NC_NOWRITE, &_ncid) );
    if (this->_statusHasError()) return;

    // dimension ids
    _status.push_back( nc_inq_dimid(_ncid, "number_of_volumes", &_nvolumesid) );
    _status.push_back( nc_inq_dimid(_ncid, "number_of_vertices", &_nverticesid) );
    _status.push_back( nc_inq_dimid(_ncid, "number_of_points", &_npointsid) );
    _status.push_back( nc_inq_dimid(_ncid, "number_of_timesteps", &_ntimestepsid) );
    if (this->_statusHasError()) return;

    // dimension values
    _status.push_back( nc_inq_dimlen(_ncid, _nvolumesid, &_nvolumes) );
    _status.push_back( nc_inq_dimlen(_ncid, _nverticesid, &_nvertices) );
    _status.push_back( nc_inq_dimlen(_ncid, _npointsid, &_npoints) );
    _status.push_back( nc_inq_dimlen(_ncid, _ntimestepsid, &_ntimesteps) );
    if (this->_statusHasError()) return;

    // debug info
    osg::notify(osg::INFO) << "[SWWReader]" <<  std::endl;
    osg::notify(osg::INFO) << "    nvolumes: " << _nvolumes <<  std::endl;
    osg::notify(osg::INFO) << "    nvertices: " << _nvertices << std::endl;
    osg::notify(osg::INFO) << "    npoints: " << _npoints << std::endl;
    osg::notify(osg::INFO) << "    ntimesteps: " << _ntimesteps << std::endl;


    // variable ids
    _status.push_back( nc_inq_varid (_ncid, "x", &_xid) );
    _status.push_back( nc_inq_varid (_ncid, "y", &_yid) );
    _status.push_back( nc_inq_varid (_ncid, "z", &_zid) );
    _status.push_back( nc_inq_varid (_ncid, "volumes", &_volumesid) );
    _status.push_back( nc_inq_varid (_ncid, "time", &_timeid) );
    _status.push_back( nc_inq_varid (_ncid, "stage", &_stageid) );
    if (this->_statusHasError()) return;

    // allocation of variable arrays, destructor responsible for cleanup
    _px = new float[_npoints];
    _py = new float[_npoints];
    _pz = new float[_npoints];
    _ptime = new float[_ntimesteps];
    _pvolumes = new unsigned int[_nvertices * _nvolumes];
    _pstage = new float[_npoints];

    // loading variables from netcdf file
    _status.push_back( nc_get_var_float (_ncid, _xid, _px) );  // x vertices
    _status.push_back( nc_get_var_float (_ncid, _yid, _py) );  // y vertices
    _status.push_back( nc_get_var_float (_ncid, _zid, _pz) );  // bedslope heights
    _status.push_back( nc_get_var_float (_ncid, _timeid, _ptime) );  // time array
    _status.push_back( nc_get_var_int (_ncid, _volumesid, (int *) _pvolumes) );  // triangle indices
    if (this->_statusHasError()) return;

    // loop index
    size_t iv;

    // vertex indices
    unsigned int v1index, v2index, v3index;

    // culling defaults
    _cullnearzero = true;
    _cullnearzerovalue = 0.001;
    _cullsteepangle = true;
    _cullsteepanglevalue = 0.9;

    // compute triangle connectivity, a list (indexed by vertex number) 
    // of lists (indices of triangles sharing this vertex)
    _connectivity = std::vector<triangle_list>(_npoints);
    for (iv=0; iv < _nvolumes; iv++)
    {
        v1index = _pvolumes[3*iv+0];
        v2index = _pvolumes[3*iv+1];
        v3index = _pvolumes[3*iv+2];
        _connectivity.at(v1index).push_back(iv);
        _connectivity.at(v2index).push_back(iv);
        _connectivity.at(v3index).push_back(iv);
    }

    // load bedslope vertex array
    _bedslopevertices = new osg::Vec3Array;
    for (iv=0; iv < _npoints; iv++)
        _bedslopevertices->push_back( osg::Vec3(_px[iv], _py[iv], _pz[iv]) );

    // load bedslope index array, pvolumes array indexes into x, y and z
    _bedslopeindices = new osg::UIntArray;
    for (iv=0; iv < _nvolumes*_nvertices; iv++)
        _bedslopeindices->push_back( _pvolumes[iv] );

    // calculate bedslope primitive normal and centroid arrays
    _bedslopenormals = new osg::Vec3Array;
    _bedslopecentroids = new osg::Vec3Array;
    osg::Vec3 v1, v2, v3, side1, side2, nrm;
    for (iv=0; iv < _nvolumes; iv++)
    {
        v1index = _pvolumes[3*iv+0];
        v2index = _pvolumes[3*iv+1];
        v3index = _pvolumes[3*iv+2];

        v1 = _bedslopevertices->at(v1index);
        v2 = _bedslopevertices->at(v2index);
        v3 = _bedslopevertices->at(v3index);

        side1 = v2 - v1;
        side2 = v3 - v2;
        nrm = side1^side2;
        nrm.normalize();

        _bedslopenormals->push_back( nrm );
        _bedslopecentroids->push_back( (v1+v2+v3)/3.0 );
    }

    // success!
    _valid = true;

}



SWWReader::~SWWReader()
{
    _status.push_back( nc_close(_ncid) );
    delete[]( _px, _py, _pz, _ptime, _pvolumes, _pstage );
}



osg::ref_ptr<osg::Vec3Array> SWWReader::getStageVertexArray(unsigned int index)
{
    size_t start[2], count[2], iv;
    const ptrdiff_t stride[2] = {1,1};
    start[0] = index;
    start[1] = 0;
    count[0] = 1;
    count[1] = _npoints;

    // stage heights from netcdf file (x and y are same as bedslope)
    int status = nc_get_vars_float (_ncid, _stageid, start, count, stride, _pstage);

    // load stage vertex array
    _stagevertices = new osg::Vec3Array;
    for (iv=0; iv < _npoints; iv++)
        _stagevertices->push_back( osg::Vec3(_px[iv], _py[iv], _pstage[iv]) );

    // stage index and per primitive normal and centroid arrays
    _stageprimitivenormals = new osg::Vec3Array;
    _stagecentroids = new osg::Vec3Array;
    _stageindices = new osg::UIntArray;
    osg::Vec3 v1b, v2b, v3b;
    osg::Vec3 v1s, v2s, v3s;
    osg::Vec3 side1, side2, nrm;
    unsigned int v1index, v2index, v3index;
    osg::ref_ptr<osg::UIntArray> culledlist = new osg::UIntArray;
    std::vector<triangle_list> stageconnectivity = std::vector<triangle_list>(_npoints);

    // over all stage triangles
    unsigned int nonculled = 0;
    for (iv=0; iv < _nvolumes; iv++){

        v1index = _pvolumes[3*iv+0];
        v2index = _pvolumes[3*iv+1];
        v3index = _pvolumes[3*iv+2];

        v1s = _stagevertices->at(v1index);
        v2s = _stagevertices->at(v2index);
        v3s = _stagevertices->at(v3index);
        v1b = _bedslopevertices->at(v1index);
        v2b = _bedslopevertices->at(v2index);
        v3b = _bedslopevertices->at(v3index);

        // shallow water depth culling test
        if( _cullnearzero && ((v1s.z()+v2s.z()+v3s.z())/3.0 - _bedslopecentroids->at(iv).z()) < _cullnearzerovalue )
        {
            culledlist->push_back( iv );
            continue;
        }

        // steep bedslope culling test
        if( _cullsteepangle )
        {
            if( _bedslopenormals->at(iv) * osg::Vec3f(0,0,1) < _cullsteepanglevalue )
            {
                culledlist->push_back( iv );
                continue;
            }
        }

        // pass through, current triangle is visible (not culled)
        _stageindices->push_back( v1index );
        _stageindices->push_back( v2index );
        _stageindices->push_back( v3index );

        // current triangle primitive normal
        side1 = v2s - v1s;
        side2 = v3s - v2s;
        nrm = side1^side2;
        nrm.normalize();

        // stage triangle connectivity, a list (indexed by vertex number) 
        // of lists (indices of non-culled triangles sharing this vertex)
        stageconnectivity.at(v1index).push_back(nonculled);
        stageconnectivity.at(v2index).push_back(nonculled);
        stageconnectivity.at(v3index).push_back(nonculled);

        // store centroid and primitive normal (both needed to visualize vectors)
        _stagecentroids->push_back( (v1s+v2s+v3s)/3.0 );
        _stageprimitivenormals->push_back( nrm );

        // non-culled triangle count
        nonculled++;

    }

    // per-vertex normals
    _stagevertexnormals = new osg::Vec3Array;
    int num_shared_triangles, triangle_index;
    for (iv=0; iv < _npoints; iv++)
    {
        nrm.set(0,0,0);
        num_shared_triangles = stageconnectivity.at(iv).size();
        for (int i=0; i<num_shared_triangles; i++ )
        {
            triangle_index = stageconnectivity.at(iv).at(i);

            // summing of contributing primitive normals
            nrm += _stageprimitivenormals->at(triangle_index);
        }

        // FIXME: Length of this list has to match vertex list length.
        // We have vertices and normals that aren't being referenced.
        if( 1 || num_shared_triangles )  // avoid vertices belonging only to culled triangles
        {
            nrm = nrm / num_shared_triangles;  // average
            nrm.normalize();
            _stagevertexnormals->push_back(nrm);
        }
    }

    return _stagevertices;
}


bool SWWReader::_statusHasError()
{
    bool haserror = false;  // assume success, trap failure

    std::vector<int>::iterator iter;
    for (iter=_status.begin(); iter != _status.end(); iter++)
    {
        if (*iter != NC_NOERR)
        {
            osg::notify(osg::WARN) << nc_strerror(*iter) <<  std::endl;
            haserror = true;
            nc_close(_ncid);
            break;
        }
    }

    // on return start gathering result values afresh
    _status.clear();

    return haserror;
}