igbhead.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/>.
// ----------------------------------------------------------------------------


/**********************************************************************

  igbhead - program to display and remove/modify/create IGB headers

  usage igbhead [options] file

  author: Edward Vigmond

 **********************************************************************/
#include<cstdio>
#include<string>
#include<iostream>
#include<libgen.h>
#include "IGBheader.h"
#include "cmdline.h"

using namespace std;

void output_header( IGBheader* );
int  jive(IGBheader *header, gzFile in );
bool iswrite(int i );
void write_slice(IGBheader *header, IGBheader *in, string ofnb, double *buf, int);
bool iscompressed (string fn);

class FILE_PTR {
  public:
    FILE_PTR( gzFile g ):gzip(true),gp(g){}
    FILE_PTR( FILE *f ):gzip(false),fp(f){}
    void close(){ if ( gzip )gzclose(gp); else fclose(fp); }
    void write(void *buf, int n){if(gzip) gzwrite( gp, buf, n );
                                 else fwrite( buf, n, 1, fp ); }
    void printfl(float f){if(gzip)gzprintf(gp,"%.7f\n",f);else
                      fprintf( fp, "%.7f\n", f ); }
  private:
        bool    gzip;
    FILE   *fp;
    gzFile  gp;
};


void write_to_binary( IGBheader *h, gzFile in, FILE_PTR *out_fp, int skip )
{
   char *buffer = new char[h->slice_sz()*h->data_size()];

   int  nr;
   int  tm = 0;
   char lbuf[2048];

   for( int i=0; i<skip; i++ )
     gzgets( in, lbuf, 2048 ); 

   do {
     nr = 0;
     double d[9];
     for( int i=0; i<h->slice_sz(); i++ ) {
       if( !gzgets( in, lbuf, 2048 ) )
         break;
       if( sscanf( lbuf, "%lf %lf %lf %lf %lf %lf %lf %lf %lf", d, d+1,
             d+2, d+3, d+4, d+5, d+6, d+7, d+8 ) != h->num_components() )
           break;
       nr++;
       h->to_bin( buffer+i*h->data_size(), d );
     }
     if( nr==h->slice_sz() )
       out_fp->write( buffer,h->slice_sz()*h->data_size() );
   }while( nr==h->slice_sz() && ++tm<h->t() );

   if( tm != h->t() )
     fprintf( stderr, "\nTimes do not match!!!!!!!!!!!\n\n" );
  
  delete[] buffer;
}


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

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

  if( argc == 1 || args_info.inputs_num != 1 ) {
    cmdline_parser_print_help();
  exit(0);
  }

  // the last arg must be a file name
  gzFile in = gzopen( args_info.inputs[0], "r" );
  if( in == NULL ) {
    cerr << "File not found: " << args_info.inputs[0] << endl;
    exit(1);
  }
  IGBheader* head_in= new IGBheader( in );
  if( head_in->read(true) != 0 ) gzrewind( in );

  // just output the current header
  if( argc==2 ) {
  output_header( head_in );
  gzclose( in );
  exit(0);
  }

  IGBheader* head_out = args_info.convert_data_flag ? new IGBheader(head_in) : head_in;

  if( args_info.x_given ) head_out->x( static_cast<size_t>(args_info.x_arg) );
  if( args_info.y_given ) head_out->y( static_cast<size_t>(args_info.y_arg) );
  if( args_info.z_given ) head_out->z( static_cast<size_t>(args_info.z_arg) );
  if( args_info.t_given ) head_out->t( static_cast<size_t>(args_info.t_arg) );
  if( args_info.data_type_given ) head_out->type( args_info.data_type_arg );
  if( args_info.system_given ) head_out->systeme( args_info.system_arg==system_arg_big ? "big_endian" : "little_endian" );
  if( args_info.dim_x_given ) head_out->dim_x( args_info.dim_x_arg );
  if( args_info.dim_y_given ) head_out->dim_y( args_info.dim_y_arg );
  if( args_info.dim_z_given ) head_out->dim_z( args_info.dim_z_arg );
  if( args_info.dim_t_given ) head_out->dim_t( args_info.dim_t_arg );
  if( args_info.org_x_given ) head_out->org_x( args_info.org_x_arg );
  if( args_info.org_y_given ) head_out->org_y( args_info.org_y_arg );
  if( args_info.org_z_given ) head_out->org_z( args_info.org_z_arg );
  if( args_info.org_t_given ) head_out->org_t( args_info.org_t_arg );
  if( args_info.inc_x_given ) head_out->inc_x( args_info.inc_x_arg );
  if( args_info.inc_y_given ) head_out->inc_y( args_info.inc_y_arg );
  if( args_info.inc_z_given ) head_out->inc_z( args_info.inc_z_arg );
  if( args_info.inc_t_given ) head_out->inc_t( args_info.inc_t_arg );
  if( args_info.x_units_given ) head_out->unites_x( args_info.x_units_arg );
  if( args_info.y_units_given ) head_out->unites_y( args_info.y_units_arg );
  if( args_info.z_units_given ) head_out->unites_z( args_info.z_units_arg );
  if( args_info.t_units_given ) head_out->unites_t( args_info.t_units_arg );
  if( args_info.clear_comment_given ) head_out->comment(NULL);
  if( args_info.comment_given ) head_out->comment( args_info.comment_arg );
  if( args_info.data_factor_given ) head_out->facteur(args_info.data_factor_arg);
  if( args_info.data_zero_given ) head_out->zero( args_info.data_zero_arg );
  if( args_info.author_given ) head_out->aut_name( args_info.author_arg );
  if( args_info.transparent_given ){
    if( strlen(args_info.transparent_arg)!=2*Data_Size[head_out->type()] ){
      cerr << "Incorrect tranparent value specified\n";
      exit(1);
    }
    // convert hex digits to bytes
    char* v = new char[Data_Size[head_out->type()]];
    char s[3], *pp;
    s[2] = '\0';
    for( int i=0; i<Data_Size[head_out->type()]; i++ ){
      s[0] = args_info.transparent_arg[i*2];
      s[1] = args_info.transparent_arg[i*2+1];
      v[i] = strtol( s, &pp, 16 );
    }
    head_out->transparent( v );
  }
  if( args_info.no_transparent_given ) head_out->transparent( NULL );

  // count the number of data bytes
  if( args_info.jive_time_given )
    jive( head_in, in );

  // make a temporary file
  string tmpfn = ".";
  tmpfn += basename(args_info.inputs[0]);
  tmpfn += ".tmp";

  // determine output file
  string ofn;
  if( args_info.frankenstein_given ){
    ofn = args_info.frankenstein_arg;
    gzclose( in );
    IGBheader *h = new IGBheader(in=gzopen( args_info.frankenstein_arg,"r" ));
    if( h->read() != 0 ) gzrewind( in );
    delete h;
  }
  else
    ofn = args_info.inputs[0];

  if(args_info.output_file_given)
    ofn = args_info.output_file_arg;

  FILE_PTR *out_fp;
  if( iscompressed(ofn) ) {
    gzFile gzout;
  head_out->fileptr(gzout = gzopen( tmpfn.c_str(), "w" ) );
    out_fp = new FILE_PTR( gzout );
  } else {
    FILE* fout;
  head_out->fileptr(fout= fopen( tmpfn.c_str(), "w" ) );
    out_fp = new FILE_PTR( fout );
  }

  if( !args_info.decapitate_given )
    head_out->write();

  if( args_info.create_given )
    write_to_binary( head_out, in, out_fp, args_info.create_arg );
  else if( args_info.convert_data_flag ) {
    double *slice_data = new double[head_in->slice_sz()];
    for( int t=0; t<head_in->t(); t++ ) {
      head_in->read_data( slice_data );
      head_out->write_data( slice_data );
    }
    delete [] slice_data;
  } else if( args_info.transpose_flag ) {
    int ds        = head_out->data_size();
    long offset   = head_out->t()*ds;
    long datStart = gztell( in );
    char *buf = (char *)malloc( head_out->slice_sz()*ds );
    for( int i=0; i<head_out->t(); i++ ) {
      for( int j=0; j<head_out->slice_sz(); j++ ) {
        gzseek( in,datStart+offset*j+i*ds, SEEK_SET );
        gzread( in, buf+j*ds, ds*1  );
      }
      out_fp->write( buf, ds*head_out->slice_sz() );
    }
    delete [] buf;
  } else {
    const int bufsize=1000000;
    char* buf = new char[bufsize];
    int nb;
    while( (nb=gzread( in, buf, bufsize )) > 0 ) out_fp->write( buf, nb );
    delete [] buf;
  }

  out_fp->close();
  gzclose( in );
  rename(  tmpfn.c_str(), ofn.c_str() );
  delete out_fp;
  return 0;
}

void output_header( IGBheader* header )
{
  bool tf;
    printf( "x dimension:\t%zu\n", header->x() );
    printf( "y dimension:\t%zu\n", header->y() );
    printf( "z dimension:\t%zu\n", header->z() );
    printf( "t dimension:\t%zu\n", header->t() );
    printf( "data type:\t%s\n", Header_Type[header->type()] );
    header->unites(tf);
  if( tf ) {
        printf( "Pixel units:\t%s\n", header->unites() );
    }
  if( header->transparent() != NULL ) {
        printf( "Transparent:\t%s\n", header->transparentstr() );
  }
    header->zero(tf);
  if( tf ) {
        printf( "Pixel zero:\t%g\n", header->zero() );
    }
    header->unites(tf);
  if( tf ) {
        printf( "Pixel scaling:\t%g\n", header->facteur() );
    }
    header->dim_x(tf);
  if( tf ) {
        printf( "X size:\t\t%g\n", header->dim_x() );
    }
    header->unites_x(tf);
  if( tf ) {
        printf( "X units:\t%s\n", header->unites_x() );
    }
    header->inc_x(tf);
  if( tf ) {
        printf( "Increment in x:\t%g\n", header->inc_x() );
    }
    header->org_x(tf);
  if( tf ) {
        printf( "X origin:\t%g\n", header->org_x() );
    }
    header->dim_y(tf);
  if( tf ) {
        printf( "Y size:\t\t%g\n", header->dim_y() );
    }
    header->unites_y(tf);
  if( tf ) {
        printf( "Y units:\t%s\n", header->unites_y() );
    }
    header->inc_y(tf);
  if( tf ) {
        printf( "Increment in y:\t%g\n", header->inc_y() );
    }
    header->org_y(tf);
  if( tf ) {
        printf( "Y origin:\t%g\n", header->org_y() );
    }
    header->dim_z(tf);
  if( tf ) {
        printf( "Z size:\t\t%g\n", header->dim_z() );
    }
    header->unites_z(tf);
  if( tf ) {
        printf( "Z units:\t%s\n", header->unites_z() );
    }
    header->inc_z(tf);
  if( tf ) {
        printf( "Increment in z:\t%g\n", header->inc_z() );
    }
    header->org_z(tf);
  if( tf ) {
        printf( "Z origin:\t%g\n", header->org_z() );
    }
    header->dim_t(tf);
  if( tf ) {
        printf( "T size:\t\t%g\n", header->dim_t() );
    }
  header->unites_t(tf);
  if( tf ) {
        printf( "T units:\t%s\n", header->unites_t() );
    }
    header->inc_t(tf);
  if( tf ) {
        printf( "Increment in t:\t%g\n", header->inc_t() );
    }
    header->org_t(tf);
  if( tf ) {
        printf( "T origin:\t%g\n", header->org_t() );
    }
  printf( "Created on:\t%s\n", header->systemestr());
    header->aut_name(tf);
  if( tf ) {
        printf( "Author:\t\t%s\n", header->aut_name() );
    }
  if( header->comment() != NULL ){
      char **cm = header->comment();
    int i=0;
    while( cm[i] != NULL )
          printf( "#%s\n", cm[i++] );
    }
}

int jive(IGBheader *header, gzFile in )
{
  z_off_t zo = gztell( in );
  const int bufsize=8196;
  char buff[bufsize];
  long nb=0, nr;
  while( (nr=gzread( in, buff, bufsize )) == bufsize )
    nb += nr;
  nb += nr;
  int slicesize = header->x()*header->y()*header->z()*header->data_size();
  header->t( nb/slicesize );
  gzseek( in, zo, SEEK_SET );
  if( header->dim_t() )
    header->dim_t( (header->t()-1)*header->inc_t() );
  return(slicesize);
}

void write_slice( IGBheader *hout, IGBheader *hin, string ofnb, double *buf, int slnum )
{
  long   nr;

  string ofn = ofnb;
  if( ofn.rfind(".gz")==ofn.length()-3 ) ofn.resize(ofn.length()-3);
  if( ofn.rfind(".igb")==ofn.length()-4 ) ofn.resize(ofn.length()-3);

  char   fnm[256];
  if( (nr=hin->read_data(buf)) == hin->slice_sz() )  {

    snprintf( fnm, sizeof fnm, "%s%d.%s", ofn.c_str(), slnum, "igb" );

  if( iscompressed(ofnb) ) {
      gzFile gzout;
    hout->fileptr(gzout = gzopen( strcat(fnm,".gz"), "w" ) );
  } else {
      FILE * fout;
    hout->fileptr(fout= fopen( fnm, "w" ) );
  }

  hout->write();
  hout->write_data( buf );
  hout->close();
  } else {
    fprintf(stderr, "Not enough slices in file\n" );
    exit(-1);
  }
}

bool iscompressed (string fn)
{
  return fn.rfind(".gz") == fn.length()-3;
}