IGBfilament.cc

Go to the documentation of this file.

/*
 * Calculate phase signularities based on Martin Bishop's isovalue approach.
 */
 
#include <string>
#include "IGBheader.h"
#include "filament.h"
#include <iostream>
#include "fil_cmdline.h"
 
#define BUFSIZE  1024
 
int  initTBuffer (tbuffer *tbf, double dt);
int  shiftTBuffer(tbuffer *tbf);
void freeTBuffer (tbuffer *tbf);
void read_mesh(std::string,ElemList *, NodeList*);
 
//
#undef __FUNCT__
#define __FUNCT__ "main"
int main( int argc, char *argv[] )
{
  ElemList      el;
  ElemList     *elst = &el;
  NodeList      nl;
  NodeList     *nlst = &nl;
  char         *dn  = NULL;
  
  struct gengetopt_args_info args;
 
  if (cmdline_parser(argc, argv, &args) != 0)
      return 1;
 
  // Positional arguments are all argv after all options
    // Using cmdline_parser’s internal method:
    int first_pos = argc - args.inputs_num; // inputs_num counts leftover arguments
    if (args.inputs_num != 2) {
        fprintf(stderr, "Error: Missing required positional arguments.\n");
        cmdline_parser_print_help();
        return 1;
    }
  char *basemesh = argv[first_pos];      // first positional argument
  char *igb_data_file = argv[first_pos+1]; // second positional argument
 
  float dt  = args.dt_arg;
  float vth = args.threshold_arg;
  std::string aux_mesh = args.aux_mesh_arg;
 
 
  read_mesh(const_cast<char*>(basemesh), elst, nlst);
 
  // initialize time slice buffer
  tbuffer tbf;
  tbf.dname  = igb_data_file;
  tbf.dsize  = sizeof(float);
  tbf.rwdth  = 2;  // two samples where we resample signal
  tbf.t      = 0.;
 
  int nSliceRead = initTBuffer(&tbf,dt);
  if( nSliceRead<0 ) exit(1);
 
  // initialize filaments
  FILE     *filHdls[3];
  init_filament_IO(aux_mesh.c_str(),tbf.numSamps,filHdls);
 
  // write bogus empty frames to keep the #slices the same as the IGB file 
  for(int i=0;i<nSliceRead/2;i++)  {
    filament f;
    write_filaments(&f,filHdls,tbf.t);
    tbf.t += tbf.dt_in;
  }
 
  // loop over file
  for(int i=0;i<tbf.numSamps-nSliceRead;i++)  {
    
    log_msg(NULL,0,0, "\nAnalyzing instant %.3f ms --- ", tbf.t);
    filament  f;
    if(tbf.filled)  {
 
      // make sure data are at +/- dt/2
      interpolateTBuffer(&tbf);
 
      // dm1 and d1 hold data now at the interval dt
      SingularityDetector_t meth = face_based;
      compute_filament((float*)tbf.dr[0], (float*)tbf.dr[1], vth,elst, nlst, &f, meth);
    }
    write_filaments(&f,filHdls,tbf.t);
 
    int err = shiftTBuffer(&tbf);
  }
 
  // write bogus empty frames to keep the #slices the same as the IGB file 
  for(int i=0;i<nSliceRead-nSliceRead/2;i++)  {
    filament f;
    write_filaments(&f,filHdls,tbf.t);
    tbf.t += tbf.dt_in;
  }
  
  log_msg(NULL,0,0, "\n");
 
  // close files
  close_filament_IO(filHdls);
  fclose( tbf.data );
 
  // deconstruct
  freeTBuffer(&tbf);
 
  cmdline_parser_free(&args);
}
 
 
int
initTBuffer(tbuffer *tbf,double dt)
{
  // open data file first
  int    IOerr = 0;
  tbf->data = fopen(tbf->dname,"rb");
  if(tbf->data!=NULL)  
    tbf->igb = new IGBheader( tbf->data, true );
 
  tbf->numSamps = tbf->igb->t();
  tbf->nodes    = tbf->igb->x();
  tbf->dt_in    = tbf->igb->inc_t();
  tbf->h_wdth   = ceil(dt/tbf->dt_in/2);;
  tbf->wdth     = tbf->h_wdth*2+1;
 
  float h_range_T = tbf->h_wdth*tbf->dt_in;
  float h_dt      = dt/2;
  tbf->rdist      = (h_range_T-h_dt)/tbf->dt_in;  // rel distance to central sample
 
  // allocate tbuffer
  tbf->d = (float **)malloc(sizeof(float*)*tbf->wdth);
  for(int i=0;i<tbf->wdth;i++)
    tbf->d[i] = (float *)malloc(tbf->dsize*tbf->nodes);
 
  if(tbf->d[tbf->wdth-1]==NULL)
    log_msg(NULL,5,LOCAL, "Warning: Memory allocation error" );
 
  // allocate interpolated tbuffer
  tbf->dr = (float **)malloc(sizeof(void*)*tbf->rwdth);
  for(int i=0;i<tbf->rwdth;i++)
    tbf->dr[i] = (float *)malloc(tbf->dsize*tbf->nodes);
 
  // globalize memory allocation success or failure
  tbf->buffered = 0;
  tbf->filled   = false;
 
  // fill
  int err = 0;
  for(int i=tbf->h_wdth;i<tbf->wdth;i++)
    //if(rootReadSeq(tbf->data,tbf->dsize,tbf->nodes,(void *)(tbf->d[i]),
                   //tbf->bswp)!=tbf->nodes)  {
    if( tbf->igb->read_data( tbf->d[i], 1, NULL ) != tbf->nodes ) {
      log_msg(NULL, 5, 0, "Premature eof. Trying to read time slice %.2f", tbf->t );
      err++;
    } else  {
      tbf->buffered++;
    }
 
  return err?-1:tbf->wdth-tbf->h_wdth;
}
 
 
int
shiftTBuffer(tbuffer *tbf)
{
  int err   = 0;
 
  float *tmp = tbf->d[0];
  for(int j=0;j<tbf->wdth-1;j++) tbf->d[j] = tbf->d[j+1];
  tbf->d[tbf->wdth-1] = tmp;
 
  //if(rootReadSeq(tbf->data,tbf->dsize,tbf->nodes,(void *)tbf->d[tbf->wdth-1],
  //               tbf->bswp)!=tbf->nodes)  {
  if( tbf->igb->read_data( tbf->d[tbf->wdth-1], 1, NULL ) != tbf->nodes ) {
    err--;
    tbf->filled = false;
    tbf->buffered--;
    log_msg(NULL,3,0, "Premature eof. Trying to read time slice %.2f,", tbf->t );
  }
  else  {
    tbf->buffered++;
    if(tbf->buffered==tbf->wdth) tbf->filled = true;
  }
  tbf->t += tbf->dt_in;
 
  return err?0:1;
}
 
void
freeTBuffer(tbuffer *tbf)
{
  for(int i=0;i<tbf->rwdth;i++) free(tbf->dr[i]);
  for(int i=0;i<tbf->wdth; i++) free(tbf->d[i] );
 
  free(tbf->d );
  free(tbf->dr);
}
 
 
void 
read_mesh( std::string mesh, ElemList *elst, NodeList *nlst )
{
  FILE *ndin = fopen( (mesh+".pts").c_str(), "r" );
  FILE *elin = fopen( (mesh+".elem").c_str(), "r" );
  if( !ndin || !elin ) {
    std::cerr << "cannot read mesh: " << mesh <<"\n";
    exit(1);
  }
 
  size_t buffsz=2048;
  constexpr size_t buffsz_etype = 2048;
  char  *buff = (char *)malloc(buffsz*sizeof(char));;
  
  // read points file
  getline( &buff, &buffsz, ndin );
  sscanf( buff, "%d", &nlst->NNodes );
  nlst->Node = (Point*)malloc(nlst->NNodes*sizeof(Point));
  for( int i=0; i<nlst->NNodes; i++ ) {
    getline( &buff, &buffsz, ndin );
#ifdef REAL_IS_DOUBLE
    sscanf( buff, "%lf %lf %lf", &nlst->Node[i].x,&nlst->Node[i].y,&nlst->Node[i].z);
#else
    sscanf( buff, "%f %f %f", &nlst->Node[i].x,&nlst->Node[i].y,&nlst->Node[i].z);
#endif
  }
  fclose( ndin );
  // read elem file
  getline( &buff, &buffsz, elin );
  sscanf( buff, "%d", &elst->NElems );
  elst->Ele = (Elem*)malloc(elst->NElems*sizeof(Elem));
  size_t elign = 0;
  int j=0;
  for( int i=0; i<elst->NElems; i++ ) {
    Elem *e = elst->Ele+j;
    getline( &buff, &buffsz, elin );
    char etype[buffsz_etype];
    int p[8];
    sscanf( buff, "%s %d %d %d %d %d %d %d %d", etype, p, p+1,p+2, p+3, p+4, p+5, p+6, p+7 );
    if( !strcmp(etype,"Tr") ){
      e->type = Tri;
      e->Nnode= 3;
    } else if( !strcmp(etype,"Tt") ){
      e->type = Tetra;
      e->Nnode= 4;
    } else {
      elign++;
      continue;
    }
    e->N = (int*)malloc(e->Nnode*sizeof(int));
    memmove( e->N, p, e->Nnode*sizeof(p[0]) );
    j++;
  }
  if( elign ) {
    std::cerr << "Info: " << elign << " elements ignored\n";
    elst->NElems = j;
  }
  fclose( elin );
}