/*
    SWWReader

    Reader of SWW files from within OpenSceneGraph.

    copyright (C) 2004 Geoscience Australia
*/

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

#include <osgDB/Registry>
#include <osgDB/ReadFile>

#include <stdlib.h>
#include <stdio.h>
#include <gdal_priv.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;

    // 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;


	// sww file can optionally contain bedslope texture image filename
	size_t attlen;
	if( nc_inq_attlen(_ncid, NC_GLOBAL, "institution", &attlen) != NC_ENOTATT )
	{
		std::string filename;
		int nc_get_att_text(_ncid, NC_GLOBAL, "institution", char filename);
		setBedslopeTexture( std::string filename );
	}

    // loop index
    size_t iv;

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

    // bedslope extents and resultant scale factors
    float xmin, xmax, xrange;
    float ymin, ymax, yrange;
    float zmin, zmax, zrange;
    float aspect_ratio;
    xmin = _px[0]; xmax = _px[0];
    ymin = _py[0]; ymax = _py[0];
    zmin = _pz[0]; zmax = _pz[0];
    for( iv=1; iv < _npoints; iv++ )
    {
	if( _px[iv] < xmin ) xmin = _px[iv];
	if( _px[iv] > xmax ) xmax = _px[iv];
	if( _py[iv] < ymin ) ymin = _py[iv];
	if( _py[iv] > ymax ) ymax = _py[iv];
	if( _pz[iv] < zmin ) zmin = _pz[iv];
	if( _pz[iv] > zmax ) zmax = _pz[iv];
    }

    xrange = xmax - xmin;
    yrange = ymax - ymin;
    zrange = zmax - zmin;
    aspect_ratio = xrange/yrange;
    _xscale = 1.0 / xrange;
    _yscale = 1.0 / yrange;
    
    // handle case of a flat bed that doesn't necessarily pass through z=0
    _zscale = (zrange == 0.0) ? 1.0 : 1.0/zrange;

    _xoffset = xmin;
    _yoffset = ymin;
    _zoffset = zmin;
    _xcenter = 0.5;
    _ycenter = 0.5;
    _zcenter = 0.0;

    if( aspect_ratio > 1.0 )
    {
	_scale = 1.0 / xrange;
	_ycenter /= aspect_ratio;
    }
    else
    {
	_scale = 1.0 / yrange;
	_xcenter /= aspect_ratio;
    }

    // alpha-scaling defaults
    _alphamin = 0.8;        
    _alphamax = 1.0;
    _heightmin = 0.0;
    _heightmax = 1.0;    

    osg::notify(osg::INFO) << "[SWWReader] xmin: " << xmin <<  std::endl;
    osg::notify(osg::INFO) << "[SWWReader] xmax: " << xmax <<  std::endl;
    osg::notify(osg::INFO) << "[SWWReader] xscale: " << _xscale <<  std::endl;
    osg::notify(osg::INFO) << "[SWWReader] ymin: " << ymin <<  std::endl;
    osg::notify(osg::INFO) << "[SWWReader] ymax: " << ymax <<  std::endl;
    osg::notify(osg::INFO) << "[SWWReader] yscale: " << _yscale <<  std::endl;
    osg::notify(osg::INFO) << "[SWWReader] zmin: " << zmin <<  std::endl;
    osg::notify(osg::INFO) << "[SWWReader] zmax: " << zmax <<  std::endl;
    osg::notify(osg::INFO) << "[SWWReader] zscale: " << _zscale <<  std::endl;


    // 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);
    }

    // bedslope vertex array, shifting and scaling vertices to unit cube
    // centred about the origin
    _bedslopevertices = new osg::Vec3Array;
    for (iv=0; iv < _npoints; iv++)
	_bedslopevertices->push_back( osg::Vec3( (_px[iv]-_xoffset)*_scale - _xcenter, 
						 (_py[iv]-_yoffset)*_scale - _ycenter, 
						 (_pz[iv]-_zoffset)*_scale - _zcenter) );

    // bedslope index array, pvolumes array indexes into x, y and z
    _bedslopeindices = new osg::DrawElementsUShort(osg::PrimitiveSet::TRIANGLES, _nvolumes*_nvertices);
    for (iv=0; iv < _nvolumes*_nvertices; iv++)
	(*_bedslopeindices)[iv] = _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 );
}



void SWWReader::setBedslopeTexture( std::string filename )
{
    osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] filename: " << filename <<  std::endl;
   _bedslopetexture = new std::string(filename);
   _bedslopegeodata.hasData = false;

   // GDAL Geospatial Image Library
   GDALDataset *pGDAL;
   GDALAllRegister();
   pGDAL = (GDALDataset *) GDALOpen( filename.c_str(), GA_ReadOnly );
   if( pGDAL == NULL )
      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] GDAL unable to open file: " << filename <<  std::endl;
   else
   {
      osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] GDAL Driver: "
			     << pGDAL->GetDriver()->GetDescription() << "/" 
			     << pGDAL->GetDriver()->GetMetadataItem( GDAL_DMD_LONGNAME ) 
			     << std::endl;

      _bedslopegeodata.xresolution = pGDAL->GetRasterXSize();
      _bedslopegeodata.yresolution = pGDAL->GetRasterYSize();

      // check if image contains geodata
      double geodata[6];
      if( pGDAL->GetGeoTransform( geodata ) == CE_None )
      {
	  _bedslopegeodata.xorigin = geodata[0];
	  _bedslopegeodata.yorigin = geodata[3];
	  _bedslopegeodata.rotation = geodata[2];
	  _bedslopegeodata.xpixel = geodata[1];
	  _bedslopegeodata.ypixel = geodata[5];
	  _bedslopegeodata.hasData = true;

	  osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xorigin: " << _bedslopegeodata.xorigin <<  std::endl;
	  osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] yorigin: " << _bedslopegeodata.yorigin <<  std::endl;
	  osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] xpixel: " << _bedslopegeodata.xpixel <<  std::endl;
	  osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] ypixel: " << _bedslopegeodata.ypixel <<  std::endl;
	  osg::notify(osg::INFO) << "[SWWReader::setBedslopetexture] rotation: " << _bedslopegeodata.rotation <<  std::endl;
      }
   }

}



osg::Image* SWWReader::getBedslopeTexture()
{
   return osgDB::readImageFile( _bedslopetexture->c_str() );
}



osg::ref_ptr<osg::Vec2Array> SWWReader::getBedslopeTextureCoords()
{

   // bedslope texture coords based on (x,y) locations scaled by extents into range [0,1]
   _bedslopetexcoords = new osg::Vec2Array;
   for( unsigned int iv=0; iv < _npoints; iv++ )
      _bedslopetexcoords->push_back( osg::Vec2( (_px[iv]-_xoffset)*_xscale, (_py[iv]-_yoffset)*_yscale ) );

   return _bedslopetexcoords;
}



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, scaling and shifting vertices to lie in the unit cube
    _stagevertices = new osg::Vec3Array;
    for (iv=0; iv < _npoints; iv++)
        _stagevertices->push_back( osg::Vec3( (_px[iv]-_xoffset)*_scale - _xcenter, 
					      (_py[iv]-_yoffset)*_scale - _ycenter, 
					      (_pstage[iv]-_zoffset)*_scale - _zcenter) );

    // stage index, per primitive normal and centroid arrays
    _stageprimitivenormals = new osg::Vec3Array;
    _stagecentroids = new osg::Vec3Array;
    osg::Vec3 v1b, v2b, v3b;
    osg::Vec3 v1s, v2s, v3s;
    osg::Vec3 side1, side2, nrm;
    unsigned int v1index, v2index, v3index;

    // over all stage triangles
    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);

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

        // store centroid and primitive normal
        _stagecentroids->push_back( (v1s+v2s+v3s)/3.0 );
        _stageprimitivenormals->push_back( nrm );
    }


    // stage height above bedslope mapped as alpha value
    //      alpha = min( a(h-hmin) + alphamin, alphamax),  h >= hmin
    //      alpha = 0,                                     h < hmin
    // where a = (alphamax-alphamin)/(hmax-hmin)
    float alpha, height, alphascale;
    alphascale = (_alphamax - _alphamin) / (_heightmax - _heightmin);
    _stagecolors = new osg::Vec4Array;
    for (iv=0; iv < _npoints; iv++)
    {
	height = _stagevertices->at(iv).z() - _bedslopevertices->at(iv).z();
	
	if (height < _heightmin)
	{
	    alpha = 0.0;
	} 
	else 
	{
	    alpha = alphascale * (height - _heightmin) + _alphamin;
	    if( alpha > _alphamax ) 
		alpha = _alphamax;	  
	}
	_stagecolors->push_back( osg::Vec4( 1.0, 1.0, 1.0, alpha ) );
    }


    // per-vertex normals calculated as average of primitive normals
    // from contributing triangles
    _stagevertexnormals = new osg::Vec3Array;
    int num_shared_triangles, triangle_index;
    for (iv=0; iv < _npoints; iv++)
    {
        nrm.set(0,0,0);
        num_shared_triangles = _connectivity.at(iv).size();
        for (int i=0; i<num_shared_triangles; i++ )
        {
            triangle_index = _connectivity.at(iv).at(i);
            nrm += _stageprimitivenormals->at(triangle_index);
        }

	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;
}