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 "head_cmdline.h"
 
using namespace std;
 
void output_header( IGBheader* );
size_t jive(IGBheader *header, gzFile in );
bool iswrite(int i );
bool has_gzip_extension(const string& fn);
 
 
void write_to_binary( IGBheader *h, gzFile in, FILE *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() )
       fwrite( buffer, h->data_size(), h->slice_sz(), out_fp );
   }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;
 
  if( has_gzip_extension(ofn) ) {
    cerr << "Compressed IGB output is not supported. Use -f <name>.igb to write an uncompressed copy." << endl;
    gzclose(in);
    return 1;
  }
 
  FILE* out_fp = fopen( tmpfn.c_str(), "wb" );
  if( out_fp == NULL ) {
    cerr << "Unable to open temporary output file: " << tmpfn << endl;
    gzclose(in);
    return 1;
  }
  head_out->fileptr(out_fp);
 
  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  );
      }
      fwrite( buf, ds, head_out->slice_sz(), out_fp );
    }
    delete [] buf;
  } else {
    const int bufsize=1000000;
    char* buf = new char[bufsize];
    int nb;
    while( (nb=gzread( in, buf, bufsize )) > 0 ) fwrite( buf, 1, nb, out_fp );
    delete [] buf;
  }
 
  fclose(out_fp);
  gzclose( in );
  rename(  tmpfn.c_str(), ofn.c_str() );
  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++] );
    }
}
 
size_t 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;
  size_t 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);
}
 
bool has_gzip_extension(const string& fn)
{
  return fn.size() >= 3 && fn.compare(fn.size()-3, 3, ".gz") == 0;
}