any2matlab.cc

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <sstream>
#include <cstdio>
#include <ctype.h>
#include <cstring>
#include <string>
#include <algorithm>
#include <vector>
#include <datrwxx/readany.h>
#include <libtime++.h>
#include <aff/series.h>
#include <gsexx.h>
#include <sffxx.h>
#include "mex.h"


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{

  const char *field_names[] = {"date", "time", "tdate", "station", "channel", "auxid", "samps", "dt",
"calib", "calper", "instype", "trace"};
  int number_of_fields;
  int date_field, time_field, tdate_field, station_field, channel_field;
  int auxid_field, samps_field, dt_field, calib_field, calper_field;
  int instype_field, trace_field;
  std::vector<int> ntraces;
  ntraces.reserve(100);

  // calculate number of fields
  number_of_fields = sizeof(field_names) / sizeof(*field_names);

  // check: if input is correct 
  if (nrhs == 0)
    mexErrMsgTxt("ERROR: missing input filename!");
  // check: more than 2 arguments
  else if (nrhs > 3)
    mexErrMsgTxt("ERROR: too many input arguments!");
  // check: one argument and string
  else if ((nrhs == 1) && !mxIsChar(prhs[0])) 
    mexErrMsgTxt("ERROR: input value must be a string!");
  // check: 2 arguments and string 
  else if (nrhs == 2 && (!mxIsChar(prhs[0]) || !mxIsChar(prhs[1])))
    mexErrMsgTxt("Error: input values must be strings!");
  // input must be a row vectors 
  else if (nrhs == 2 && ((mxGetM(prhs[0]) != 1) || mxGetM(prhs[1]) != 1))
    mexErrMsgTxt("Error: input isn't valid!");
  else if (nrhs == 3 && (!mxIsChar(prhs[0]) || !mxIsChar(prhs[1]) || !mxIsChar(prhs[2])))
    mexErrMsgTxt("Error: input values must be strings!");
  else if (nrhs == 3) {
    for(int i = 0; i < 3; ++i) {
      if(mxGetM(prhs[i]) != 1) 
        mexErrMsgTxt("Error: input isn't valid!");
    }
  }

  std::string infile;
  std::string filetype;

  // If only one argument is given to any2matlab.mexa64
  if (nrhs == 1) {
    infile = std::string(mxArrayToString(prhs[0]));

  /*
        * Important note:
        * If only one argument is given to any2matlab.mexa64 the given filename must
  * contain at least one dot. After the the dot there must be one of the
  * following endings.
        */ 
        if (nrhs == 1 && infile.rfind(".") != string::npos) { 
      filetype = infile.substr(infile.rfind(".") + 1);
         }
         else  
      mexErrMsgTxt("Error: not a valid filename!");
  }

  // check: 2 arguments
  if (nrhs == 2) {
    infile = std::string(mxArrayToString(prhs[0]));
    filetype = std::string(mxArrayToString(prhs[1]));
  }

  // set format to default value: double
  std::string dtype("double");

  // check: 3 arguments
  if (nrhs == 3) {
    infile = std::string(mxArrayToString(prhs[0]));
    filetype = std::string(mxArrayToString(prhs[1]));

    std::string arg2(mxArrayToString(prhs[2]));
    // check: third argument
          if (arg2 == "int" || arg2 == "double") { 
      dtype = arg2;
          }
          else {
      mexErrMsgTxt("Error: third argument is not valid!");
    }
  }

  // set filetype to lower case
  transform(filetype.begin(), filetype.end(), filetype.begin(), tolower);
  mwSize trace = 0;
  bool hot = false;

  {      
    std::ifstream ifs(infile.c_str(), datrw::ianystream::openmode(filetype));
    datrw::ianystream is(ifs, filetype);
 
    /* get number of traces in infile */
    /* while skipping series */
    hot=is.good();

    while (hot) {
      is.skipseries();

      sff::WID2 wid2;
      // read the WID2 header
      is >> wid2;
    
      ntraces.push_back(wid2.nsamples);
      //mexPrintf("samples: %d\n", wid2.nsamples);

      trace++;

      hot=(!is.last());
    }
  }

  mwSize dims[2] = {1, trace};
  // create a 1-by-n array of strucntraces.  
  plhs[0] = mxCreateStructArray(2, dims, number_of_fields, field_names);

  date_field = mxGetFieldNumber(plhs[0], "date");
  time_field = mxGetFieldNumber(plhs[0], "time");
  tdate_field = mxGetFieldNumber(plhs[0], "tdate");
  station_field = mxGetFieldNumber(plhs[0], "station");
  auxid_field = mxGetFieldNumber(plhs[0], "auxid");
  channel_field = mxGetFieldNumber(plhs[0], "channel");
  samps_field = mxGetFieldNumber(plhs[0], "samps");
  dt_field = mxGetFieldNumber(plhs[0], "dt");
  calib_field = mxGetFieldNumber(plhs[0], "calib");
  calper_field = mxGetFieldNumber(plhs[0], "calper");
  instype_field = mxGetFieldNumber(plhs[0], "instype");
  trace_field = mxGetFieldNumber(plhs[0], "trace");

  {
    std::ifstream ifs(infile.c_str(), datrw::ianystream::openmode(filetype));
    datrw::ianystream is(ifs, filetype);
    hot = is.good();

    int i = 0;

    while(hot) {
      mxArray *field_value;
      
      // read integer data
      if (dtype == "int") {
        datrw::Tiseries iseries;

        // read the series
        is >> iseries;

        // set trace-field
        field_value = mxCreateNumericMatrix(ntraces[i], 1, mxINT32_CLASS,
          mxREAL);
        // copy stuff
        int *p = iseries.pointer();
        int *ptr = static_cast<int*>(mxGetData(field_value));
        for(int j = 0; j < iseries.size(); ++j) {
          *ptr = *p;
          p++;
          ptr++;
        }
        mxSetFieldByNumber(plhs[0], i, trace_field, field_value);
      }
      // read double data - default
      else {
        datrw::Tdseries dseries;
        // read the series
        is >> dseries;

        // set trace-field
        field_value = mxCreateNumericMatrix(ntraces[i], 1, mxDOUBLE_CLASS,
          mxREAL);
        // copy stuff
        double *ptr = mxGetPr(field_value); 
        double *tmp = dseries.pointer();
        for(int j = 0; j < dseries.size(); ++j) {
          ptr[j] = *tmp;
          tmp++;
        }
        mxSetFieldByNumber(plhs[0], i, trace_field, field_value);
      }

      sff::WID2 wid2;
      // read the WID2 header
      is >> wid2;

      // set date-field
      {
        std::ostringstream str;
        str << wid2.date.year() << "/" 
          << wid2.date.month() << "/" 
          << wid2.date.day() << std::ends;
        std::string datestring; 
        datestring.append(str.str());
        mxSetFieldByNumber(plhs[0], i, date_field,
          mxCreateString(datestring.c_str()));
      }

      // set date-field
      {
        std::ostringstream str;
        std::string timestring; 
        str << wid2.date.hour() << ":" 
          << wid2.date.minute() << ":" 
          << wid2.date.second() << "." 
          << wid2.date.milsec() << std::ends;
        timestring.append(str.str());
        mxSetFieldByNumber(plhs[0], i, time_field,
          mxCreateString(timestring.c_str()));
      }

      // set tdate-field
      double time [6];
      time [0] = wid2.date.year();
      time [1] = wid2.date.month();
      time [2] = wid2.date.day();
      time [3] = wid2.date.hour();
      time [4] = wid2.date.minute();
      time [5] = wid2.date.second() + wid2.date.milsec() * 0.001;

      field_value = mxCreateDoubleMatrix(1,
        sizeof(time) / sizeof(time[0]), mxREAL);
      memmove(mxGetPr(field_value), &time,
        sizeof(time) / sizeof(time[0]) * sizeof(double));
      mxSetFieldByNumber(plhs[0], i, tdate_field, field_value);

      // set station-field
      mxSetFieldByNumber(plhs[0], i, station_field,
        mxCreateString(wid2.station.c_str()));

      // set channel-field
      mxSetFieldByNumber(plhs[0], i, channel_field,
        mxCreateString(wid2.channel.c_str()));

      // set auxid-field
      mxSetFieldByNumber(plhs[0], i, auxid_field,
        mxCreateString(wid2.auxid.c_str()));

      // set samps-field
      field_value = mxCreateDoubleMatrix(1, 1, mxREAL);
      *mxGetPr(field_value) = wid2.nsamples;
      mxSetFieldByNumber(plhs[0], i, samps_field, field_value);

      // set dt-field
      field_value = mxCreateDoubleMatrix(1, 1, mxREAL);
      *mxGetPr(field_value) = wid2.dt;
      mxSetFieldByNumber(plhs[0], i, dt_field, field_value);

      // set calib-field
      field_value = mxCreateDoubleMatrix(1, 1, mxREAL);
      *mxGetPr(field_value) = wid2.calib;
      mxSetFieldByNumber(plhs[0], i, calib_field, field_value);

      // set calper-field
      field_value = mxCreateDoubleMatrix(1, 1, mxREAL);
      *mxGetPr(field_value) = wid2.calper;
      mxSetFieldByNumber(plhs[0], i, calper_field, field_value);

      // set instype-field
      mxSetFieldByNumber(plhs[0], i, instype_field,
        mxCreateString(wid2.instype.c_str()));

      i++;
      hot=(!is.last());
    }
  }
}