IGBextract.cc

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// 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 3 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, see <https://www.gnu.org/licenses/>.
// ----------------------------------------------------------------------------


#include<stdlib.h>
#include<stdio.h>
#include<iostream>
#include "IGBheader.h"
#include <assert.h>
#include "ext_cmdline.h"
#include <vector>
#include <algorithm>

using namespace std;

bool
process_specifier( char *sp, vector<int> &nodes, int max )
{
  int inc = 1;
  if( strchr( sp, ':' ) ){
    char *colptr = strchr( sp, ':' );
    int consumed;
    sscanf( colptr+1, "%d%n", &inc, &consumed );
    if( !consumed || consumed!=strlen(colptr+1) ) {
      fprintf( stderr, "Illegal range specified: %s\n", sp );
      exit(1);
    }
    *colptr = '\0';
  }

  if( !strcmp( "*", sp ) ) {
    nodes.clear();
    for( int i=0; i<max; i+=inc )
      nodes.push_back(i);
    return inc==1;
  }

  if( strchr( sp, '-' ) ){
    int  s;
    char fs[100];
    int  nr = sscanf( sp, "%d-%s", &s, fs );
    if( nr<2 ) {
      fprintf( stderr, "Illegal range specified: %s\n", sp );
      exit(1);
    }

    int f;
    if( !strcmp(fs, "*") )
      f = max-1;
    else {
      int consumed;
      sscanf( fs, "%d%n", &f, &consumed );
      if( consumed != strlen(fs) ) {
        fprintf( stderr, "Illegal range specified: %s\n", sp );
        exit(1);
      }
    }

    if( s<0 || f<s ) {
      fprintf( stderr, "Illegal range specified: %s\n", sp );
      exit(1);
    }
    for( int i=s; i<=f; i+=inc )
      nodes.push_back( i );

  } else {

    int num;
    sscanf( sp, "%d", &num );

    if( !sscanf( sp, "%d", &num ) || num<0 ) {
      fprintf( stderr, "Illegal node number specified in list\n" );
      exit(1);
    }
    nodes.push_back( num );
  }
  return false;
}


void
extract_nodes( const char *list, vector<int> &nodes, int max )
{
  char *lc = strdup( list );
  char *sp = strtok( lc, "," );
  while( sp != NULL ) {

    if( process_specifier( sp, nodes, max ) ) return;

    sp = strtok( NULL, "," );
  }
  free( lc );
}


void
read_nodes( char *file, vector<int> &nodes, int max )
{
  FILE *in = fopen( file, "r" );
  if( in == NULL ) {
      fprintf( stderr, "cannot open node list file\n" );
      exit(1);
  }
#define BUFSIZE 1024
  char buf[BUFSIZE];
  while( fgets( buf, BUFSIZE, in ) != NULL ) {
      if( process_specifier( buf, nodes, max ) )
        return;
  }
  fclose( in );
}


string
outfile( string specified, enum_format format, int s=-1 )
{
  if( specified=="-" )
    return specified.c_str();

  string ext;
  switch( format ) {
      case format_arg_binary:
      case format_arg_asciiTm:
      case format_arg_ascii:
      case format_arg_ascii_1cpL:
      case format_arg_ascii_1pL:
      case format_arg_ascii_1pLn:
        // no standard extension
        break;
      case format_arg_dat_t:
        ext = ".dat_t";
        break;
      case format_arg_IGB:
        ext = ".igb";
        break;
  }

  string slice;
  if( s>=0 ) {
    char paddedS[10];
    sprintf( paddedS, "_%05d", s );
    slice = paddedS;
  }

  // remove possible extension
  if( specified.length()>ext.length() &&
          specified.substr(specified.length()-ext.length()) == ext )
    specified = specified.substr(0,specified.length()-ext.length());

  return specified + slice + ext;
}


void
output( vector<int> nodes, IGBheader *h,  char *ofname,
                 enum_format of, int t0, int t1, int stride, bool explode, float scale)
{
  IGBheader hout = *h;
  FILE     *out=NULL;
  if( !strcmp(ofname, "-") )
    out = stdout;
  else if( !explode ) {
    string out_fname = outfile(ofname, of);
    out = fopen(out_fname.c_str(), "w");
    if (!out) {
      fprintf(stderr, "Can't open output file '%s'\n", out_fname.c_str());
      exit(1);
    }
  }

  if (t1 >= h->t()) t1 = h->t()-1;
  int tend = t1==-1? h->t()-1 : t1;

  hout.x( nodes.size() );
  hout.y(1);
  hout.z(1);
  hout.t( explode ? 1 : (tend-t0)/stride+1 );
  hout.org_t( h->org_t()+t0*h->inc_t() );
  hout.inc_t(h->inc_t()*stride);
  hout.dim_t( (t1-t0)*hout.inc_t() );
  if( h->num_components()==3 )
    hout.type( IGB_VEC3_f );
  else if( h->num_components()==4 )
    hout.type( IGB_VEC4_f );
  else if( h->num_components()==9 )
    hout.type( IGB_VEC9_f );
  else
    hout.type(IGB_FLOAT);

  if( of==format_arg_IGB ) {
    hout.fileptr( out );
    if( !explode ) hout.write();
  }
  const bool lessnodes = nodes.size() !=  (size_t) h->x();
  std::vector<float> slice(h->slice_sz() * hout.num_components());
  std::vector<float> slice_cropped(lessnodes ? nodes.size() * hout.num_components() : 0);
  float * slice_ptr = slice.data();
  float *output_slice = lessnodes ? slice_cropped.data() : slice_ptr;
  
  if( of==format_arg_dat_t )
    fprintf( out, "%d\n", (tend-t0)/stride+1 );
  else if( of==format_arg_ascii_1pLn )
    fprintf( out, "%d\n", int(nodes.size()) );

  for( int t=t0; t<=tend; t+=stride ) {

    if( explode && out!=stdout ) {
      string out_fname = outfile(ofname, of, t);
      out = fopen(out_fname.c_str(), "w");
      if( of==format_arg_IGB ) {
        hout.fileptr(out);
        hout.write();
      }
    }

    if( of==format_arg_dat_t )
      fprintf(out, "%d\n", int(nodes.size()));

    h->go2slice(t);
    h->read_data( slice_ptr );

    if( of== format_arg_asciiTm )
      fprintf( out, "%f", h->org_t()+h->inc_t()*t );

    if(lessnodes) {
      size_t cnt = 0; 
      for(const auto p : nodes)
        for(size_t k=0; k < hout.num_components(); k++)
          slice_cropped[cnt++] = slice[hout.num_components() * p + k];
    }
    std::for_each(output_slice, output_slice + nodes.size() * hout.num_components(), [scale](float & c){ c *= scale; });
    switch(of) {
      case format_arg_IGB:
      case format_arg_binary:
        fwrite(output_slice, sizeof(float), hout.slice_sz() * hout.num_components(), out );
        break;
      case format_arg_asciiTm:
      {
        size_t cnt=0;
        for(size_t i = 0; i < nodes.size(); i++)
          for(size_t k=0; k < hout.num_components(); k++)
            fprintf(out, " %g", output_slice[cnt++]);
        break;
      }
      case format_arg_ascii_1cpL:
      {
        size_t cnt=0;
        for(size_t i = 0; i < nodes.size(); i++)
          for(size_t k=0; k < hout.num_components(); k++)
            fprintf(out, " %g\n", output_slice[cnt++]);
        break;
      }
      case format_arg_ascii_1pL:
      case format_arg_ascii_1pLn:
      {
        size_t cnt=0;
        for(size_t i = 0; i < nodes.size(); i++) {
          for(size_t k=0; k < hout.num_components(); k++)
            fprintf(out, " %g", output_slice[cnt++]);
          fputs("\n", out);
        }
        break;
      }
      case format_arg_ascii:
      {
        size_t cnt=0;
        for(size_t i=0; i < nodes.size(); i++) 
          for(size_t k=0; k < hout.num_components(); k++)
            fprintf(out, "%g ", output_slice[cnt++]);  
        break;
      }
      case format_arg_dat_t:
      {
        size_t cnt=0;
        for(size_t i=0; i < nodes.size(); i++)
          for(size_t k=0; k < hout.num_components(); k++)
            fprintf(out, "%g\n", output_slice[cnt++]);
        break;
      }
    }
    if( of==format_arg_ascii || of==format_arg_asciiTm )
      fprintf( out, "\n" );
    if( explode && out!=stdout )
      fclose(out);
  }

  if( !explode && out != stdout )
  fclose( out );
}


int
main( int argc, char *argv[] )
{
  gengetopt_args_info args;

  // let's call our cmdline parser
  if (cmdline_parser (argc, argv, &args) != 0)
     exit(1);

  // the last arg must be a file name
  gzFile in = ( !args.inputs_num || !strcmp(args.inputs[0],"-") ) ?
               gzdopen( dup(0), "r" ):       // duplicate stdin
               gzopen( args.inputs[0], "r" );
  if( in == NULL ) {
    cerr << "File not found: " << args.inputs[0] << endl;
    exit(1);
  }
  IGBheader header(in);
  header.read(true);

  if( args.explode_flag && args.format_arg==format_arg_dat_t ) {
    cerr << "Cannot explode into dat_t format" << endl;
    exit(1);
  }

  vector<int> nodes;
  if( args.list_given )
    extract_nodes( args.list_arg, nodes,  header.slice_sz() );
  else if( args.input_file_given )
    read_nodes( args.input_file_arg, nodes, header.slice_sz() );
  else // read them all
    extract_nodes( "*", nodes, header.slice_sz() );

  for( auto  p : nodes )
    if( header.slice_sz() <= p ) {
    cerr << "Error: out of range nodes requested" << endl;
    exit(1);
  }

  int f0=args.f0_arg, f1=args.f1_arg;
  if( args.t0_given ) 
    f0 = std::round((args.t0_arg-header.org_t())/header.inc_t());
  if( f0<0 || f0>=header.t() ) {
      cerr << "Illegal ending frame specified: "<<f0<<" not in [0,"<<header.t()<<")"<< endl;
      exit(1);
  }
  if( args.t1_given )
    f1 = args.t1_arg==-1 ? header.t()-1 : std::round((args.t1_arg-header.org_t())/header.inc_t());
  if( f1>=0 && f1<f0 ) {
      cerr << "Illegal ending frame specified: "<<f1<<" not in ["<<f0<<","<<header.t()<<")"<< endl;
      exit(1);
  }

  output(nodes, &header, args.output_file_arg,args.format_arg,
                                   f0, f1, args.df_arg, args.explode_flag, args.scale_arg );

  cmdline_parser_free(&args);
}