/home/pascal/depot/filedata/src/hdffiledata.h

/***************************************************************************
 *   Copyright (C) 2005 by Nicolas PASCAL                                  *
 *   nicolas.pascal@icare.univ-lille1.fr                                   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

/***************************************************************************
 *  History :                                                              *
 *  25/07/05 : start                                                       *
 ***************************************************************************/

#ifndef HDFFILEDATA_H
#define HDFFILEDATA_H

#define HDFFILEDATA_TRACE 0

#include "filedatareader.h"
#include "Hdf_file.hpp"
#include "hdffilemetadata.h"
#include "netcdf.h"

#include <cmath>
#include <cfloat>

static const float32 VAR_FILL_VALUE=-DBL_MAX;

class HDFFileData : public FileDataReader {
    template <class T>
    T* get_value( const char * sds_name = "Height", int *start = NULL, int *stride = NULL, int *edges = NULL, int rank = -1);
    const int32 get_sds_type_code( const string & sds_name );
protected:
    Hdf_file *hdf_file;
    HDFFileMetaData *metadata;
    template <class T>
    const vector<int> get_nearest_index(const vector<T>& v_data, const T &val = T(0));
    template <class T>
    const int get_nearest_index(const vector <T> *data, const T val);
    template <class T>
    const int get_nearest_index(const T* data = NULL, const int data_size = 0, const T val = T(0));

    void init_read_write_null_input_param(const char* sds_name, int *&start, int *&stride, int *&edges, int &rank, bool * initialized_values);
    void free_read_write_allocations( const bool *are_limits_initialized, int *start, int *stride, int *edges );
    const bool check_read_write_limits(const char* sds_name, int *start, int *stride, int *edges, const int rank);
public:
    /*********************
    **   CONSTRUCTORS   **
    *********************/
    HDFFileData(const string &name="",const string &mode="r", const bool open_file = false);
    HDFFileData(const HDFFileData &hfd);
    virtual ~HDFFileData();
    HDFFileData &operator= (const HDFFileData &hfd);
    Hdf_file * get_hdf_file();
    void load_hdf_file();
    void free_hdf_file();
    const bool is_hdf_file_loaded();
    HDFFileMetaData* get_metadata();
    void load_hdf_metadata();
    void free_hdf_metadata();
    const bool is_hdf_metadata_loaded();
    virtual void get_dataset_fill_value(const string &sds_name, void* fill_value);
    vector<int> get_sds_dimension(const string &sds_name);
    virtual vector<int> get_dataset_dimension(const string &sds_name) {
        return get_sds_dimension (sds_name);
    };
    vector<int> get_vdata_dimension(const string &vdata_name);

    int get_n_dataset(){return this->get_hdf_file()->get_n_sds();};

    string get_dataset_name(int i){return this->get_hdf_file()->get_sds_name(i);};

    Hdf_sds get_dataset(string sds_name){return this->get_hdf_file()->get_sds(sds_name.c_str());};

    int get_dataset_data_type(string sds_name)
    {return this->get_hdf_file()->get_sds_data_type(sds_name.c_str());};

string get_values_attr_dataset(string sds_name,string attr_name);

    bool has_attr_dataset(string sds_name,string attr_name);

    void close_data_file(){free_hdf_file();};

    void open_data_file(){load_hdf_file();};

    string get_values_attr(string attr_name);

    string get_file_attr(string attr_name);

    bool has_attr(string attr_name);

    bool has_file_attr(string attr_name);
    /************************
    ** DATA ACCESS METHODS **
    ************************/
    template <class T>
    T get_var_value(const string sds_name, const int &y_index,const int &x_index);
    template <class T>
    T get_value_0D( const char * sds_name = "Height", int *start = NULL, int *stride = NULL, int *edges = NULL, int &rank = -1);
    template <class T>
    vector <T> *get_value_1D( const char * sds_name = "Height", int *start = NULL, int *stride = NULL, int *edges = NULL, int &rank= -1);
    template <class T>
    vector< vector <T> >* get_value_2D( const char * sds_name = "Height", int *start = NULL, int *stride = NULL, int *edges = NULL, int &rank=-1);
    virtual void * read_data(void * data, const char * sds_name, int * start = NULL, int * stride = NULL, int * edges = NULL, int rank = -1);

    template<typename DataType>
    DataType * read_data(DataType * data, const char * sds_name, int * start = NULL, int * stride = NULL, int * edges = NULL, int rank = -1);

    void* read_vdata(void* data,const char* vdata_name="", const char* vdata_field="", int start=0, int edges=-1);

    void get_fillValue(const string &sds_name, void *fillValue);

    void get_scaling(const string &sds_name, float64 &scale, float64 &offset);
};

/*************************************************/
/************ TEMPLATES DEFINITION ***************/
/*************************************************/
template <class T>
T* HDFFileData::get_value( const char * sds_name, int *start, int *stride, int *edges, int rank) {
    bool are_limits_initialized[3]; // check if start, stride or edges had been initialized in this method
    bool hdf_file_already_loaded = is_hdf_file_loaded();
    try {
        if (!hdf_file_already_loaded)
            load_hdf_file();
        // --- initialize NULL input parametres ---
        init_read_write_null_input_param(sds_name, start, stride, edges,rank, are_limits_initialized);
        // check the selection limits
        check_read_write_limits(sds_name, start,stride,edges,rank); // throw bad_index in case of problem
        // check if there is a conflict between the template type and the sds type
        if ( typeid(T).name() != Hdf_common::hdf_type_info(hdf_file->get_sds_data_type(sds_name)) ) {
            bad_type e(Hdf_common::hdf_type_info(hdf_file->get_sds_data_type(sds_name)),typeid(T).name(),sds_name);
            throw e;
        }
        // allocate a buffer to read the data
        int data_size = 1;
        for (int i = 0 ; i<rank ; ++i)
            data_size*=edges[i];
        T* hdf_data = new T[data_size];
        // read it
        hdf_file->read_sds(sds_name,hdf_data,(int32*)start,(int32*)stride,(int32*)edges);
        // free the ressources allocated by this method
        free_read_write_allocations(are_limits_initialized, start, stride, edges);
        if (!hdf_file_already_loaded)
            free_hdf_file();
        return (hdf_data);
    } catch(...) {
        // free the ressources allocated in this method
        free_read_write_allocations(are_limits_initialized , start, stride,edges);
        throw;
    }
};
template <class T>
T HDFFileData::get_value_0D( const char * sds_name, int *start, int *stride, int *edges, int &rank ) {
    T* data = get_value<T>(sds_name,start,stride,edges,rank);
    T ret=(*data);
    delete[] data;
    return ret;
}
template <class T>
vector <T> *HDFFileData::get_value_1D( const char * sds_name, int *start, int *stride, int *edges, int &rank) {
    // specific to the method
    int method_rank = 1;

    // read the values in a linear buffer
    T* data = get_value<T>(sds_name,start,stride,edges,rank);

    // initialize the return vector size
    int significant_edges_index[method_rank]; // the index of the edges > 1
    int dim_size[method_rank]; // size of the edges > 1

    int dim_size_to_init_index = 0;
    for (int i=0 ; i<rank ; ++i) {
        if  (edges[i]>1) {
            dim_size[dim_size_to_init_index]=edges[i];
            significant_edges_index[dim_size_to_init_index]=i;
            ++dim_size_to_init_index;
        }
    }
    // allocate the vector
    vector <T> *v_data = new vector <T>; // format to return : easier and faster to process
    v_data->reserve(dim_size[0]);
    // put the data in the vector : easier and faster to process
    for ( int i = 0 ; i < dim_size[0] ; ++i ) {
        v_data->push_back(data[i]);
    }

    delete[] data;
    return v_data;
}
template <class T>
vector< vector <T> >* HDFFileData::get_value_2D( const char * sds_name , int *start, int *stride, int *edges, int &rank) {
    // rank specific to the method
    int method_rank = 2;
    // read the values in a linear buffer
    T* data = get_value<T>(sds_name,start,stride,edges,rank);
    // initialize the return vector size
    int significant_edges_index[method_rank]; // the index of the edges > 1
    int dim_size[method_rank]; // size of the edges > 1
    int dim_size_to_init_index = 0;
    for (int i=0 ; i<rank ; ++i) {
        if  (edges[i]>1) {
            dim_size[dim_size_to_init_index]=edges[i];
            significant_edges_index[dim_size_to_init_index]=i;
            ++dim_size_to_init_index;
        }
    }
    // allocate the vector
    vector< vector <T> > *v_data = new vector< vector <T> >; // format to return : easier and faster to process
    v_data->reserve(dim_size[0]);
    // put the data in the vector : easier and faster to process
    for ( int i = 0 ; i < dim_size[0] ; ++i ) {
        vector <T> tmp;
        for ( int j = 0 ; j < dim_size[1] ; ++j )
            tmp.push_back(data[i*dim_size[1]+j]);
        v_data->push_back(tmp);
    }
    delete[] data;
    return v_data;
}
template <class T>
T HDFFileData::get_var_value(const string var_name, const int &y_index,const int &x_index) {
    int start[]={y_index,x_index};
    int edges[]={1,1};
    int rank=2;
    return get_value_0D<T>(var_name.c_str(),start,NULL,edges,rank);
}
template <class T>
const vector<int> HDFFileData::get_nearest_index(const vector<T> &v_data, const T &val) {
    vector<int> ret;
    if (!v_data.empty()) {
        T diff = abs((v_data[0])-val);
        T current_diff; // temporary difference -> faster
        int end = v_data.size();
        for (int i = 0; i<end ; ++i) {
            current_diff=abs((v_data[i])-val);
            if (current_diff<diff) { // found a nearest value
                diff=current_diff;
                ret.clear();
                ret.push_back(i);
            } else if (current_diff==diff) { // found one more nearest value
                ret.push_back(i);
            }
        }
    }
    return ret;
}
template <class T>
const int HDFFileData::get_nearest_index(const vector <T> *data, const T val) {
    int nearest_val_index=-1;
    T nearest_value;
    T delta;
    if (!data->empty()) {
        int length = data->size();
        nearest_val_index=0;
        nearest_value=(*data)[nearest_val_index];
        delta=abs(val-nearest_value);
        for (int i = 1 ; i<length ; ++i) {
            if ( abs(val-(*data)[i])<delta ) {
                nearest_val_index=i;
                nearest_value=(*data)[nearest_val_index];
                delta=abs(val-nearest_value);
            }
        }
    }
    return nearest_val_index;
}
template <class T>
const int HDFFileData::get_nearest_index(const T *data, const int data_size, const T val) {
    int nearest_val_index=-1;
    T nearest_value;
    T delta;
    if (data!=NULL && data_size!=0) {
        nearest_val_index=0;
        nearest_value=*(data+nearest_val_index);
        delta=abs(val-nearest_value);
        for (int i = 1 ; i<data_size ; ++i) {
            if ( abs(val-*(data+i))<delta ) {
                nearest_val_index=i;
                nearest_value=*(data+nearest_val_index);
                delta=abs(val-nearest_value);
            }
        }
    }
    return nearest_val_index;
}

// ostream& operator << (ostream& o, const HDFFileData& _hdf);

template <typename T>
T * HDFFileData::read_data(T * data, const char * sds_name, int * start, int * stride, int * edges, int rank ) {
    bool are_limits_initialized[3]; // check if start, stride or edges had been initialized in this method
    bool hdf_file_already_loaded = is_hdf_file_loaded();
    try {
        if (!hdf_file_already_loaded)
            load_hdf_file();
        // --- initialize NULL input parametres ---
        init_read_write_null_input_param(sds_name,start,stride,edges,rank,are_limits_initialized);
        // check the selection limits
        check_read_write_limits(sds_name,start,stride,edges,rank); // throw bad_index in case of problem
        // check if DataType has the same type than the sds values
        if ( typeid(T).name() != Hdf_common::hdf_type_info(hdf_file->get_sds_data_type(sds_name)) ) {
            bad_type e(Hdf_common::hdf_type_info(hdf_file->get_sds_data_type(sds_name)),typeid(T).name(),sds_name);
            throw e;
        }
        // read the desired value, in the good format
        int data_size = 1;
        for (int i = 0 ; i<rank ; ++i)
            data_size*=edges[i];
        // data allocation hasn't been done by the caller
        if (data==NULL)
            // allocate data
            data=new T[data_size];
        hdf_file->read_sds(sds_name,data,(int32*)start,(int32*)stride,(int32*)edges);
        // free the ressources allocated by this method
        free_read_write_allocations(are_limits_initialized,start,stride,edges);
        if (!hdf_file_already_loaded)
            free_hdf_file();
        return data;

    } catch(...) {
        free_read_write_allocations(are_limits_initialized,start,stride,edges);
        throw;
    }
}

#endif