openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
igbhead.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // openCARP is an open cardiac electrophysiology simulator.
3 //
4 // Copyright (C) 2020 openCARP project
5 //
6 // This program is free software: you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation, either version 3 of the License, or
9 // (at your option) any later version.
10 //
11 // This program is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with this program. If not, see <https://www.gnu.org/licenses/>.
18 // ----------------------------------------------------------------------------
19 
20 
21 /**********************************************************************
22 
23  igbhead - program to display and remove/modify/create IGB headers
24 
25  usage igbhead [options] file
26 
27  author: Edward Vigmond
28 
29  **********************************************************************/
30 #include<cstdio>
31 #include<string>
32 #include<iostream>
33 #include<libgen.h>
34 #include "IGBheader.h"
35 #include "cmdline.h"
36 
37 using namespace std;
38 
39 void output_header( IGBheader* );
40 int jive(IGBheader *header, gzFile in );
41 bool iswrite(int i );
42 void write_slice(IGBheader *header, IGBheader *in, string ofnb, double *buf, int);
43 bool iscompressed (string fn);
44 
45 class FILE_PTR {
46  public:
47  FILE_PTR( gzFile g ):gzip(true),gp(g){}
48  FILE_PTR( FILE *f ):gzip(false),fp(f){}
49  void close(){ if ( gzip )gzclose(gp); else fclose(fp); }
50  void write(void *buf, int n){if(gzip) gzwrite( gp, buf, n );
51  else fwrite( buf, n, 1, fp ); }
52  void printfl(float f){if(gzip)gzprintf(gp,"%.7f\n",f);else
53  fprintf( fp, "%.7f\n", f ); }
54  private:
55  bool gzip;
56  FILE *fp;
57  gzFile gp;
58 };
59 
60 
68 void write_to_binary( IGBheader *h, gzFile in, FILE_PTR *out_fp, int skip )
69 {
70  char *buffer = new char[h->slice_sz()*h->data_size()];
71 
72  int nr;
73  int tm = 0;
74  char lbuf[2048];
75 
76  for( int i=0; i<skip; i++ )
77  gzgets( in, lbuf, 2048 );
78 
79  do {
80  nr = 0;
81  double d[9];
82  for( int i=0; i<h->slice_sz(); i++ ) {
83  if( !gzgets( in, lbuf, 2048 ) )
84  break;
85  if( sscanf( lbuf, "%lf %lf %lf %lf %lf %lf %lf %lf %lf", d, d+1,
86  d+2, d+3, d+4, d+5, d+6, d+7, d+8 ) != h->num_components() )
87  break;
88  nr++;
89  h->to_bin( buffer+i*h->data_size(), d );
90  }
91  if( nr==h->slice_sz() )
92  out_fp->write( buffer,h->slice_sz()*h->data_size() );
93  }while( nr==h->slice_sz() && ++tm<h->t() );
94 
95  if( tm != h->t() )
96  fprintf( stderr, "\nTimes do not match!!!!!!!!!!!\n\n" );
97 
98  delete[] buffer;
99 }
100 
101 
102 int main( int argc, char* argv[] )
103 {
104  gengetopt_args_info args_info;
105 
106  // let's call our cmdline parser
107  if (cmdline_parser (argc, argv, &args_info) != 0)
108  exit(1);
109 
110  if( argc == 1 || args_info.inputs_num != 1 ) {
111  cmdline_parser_print_help();
112  exit(0);
113  }
114 
115  // the last arg must be a file name
116  gzFile in = gzopen( args_info.inputs[0], "r" );
117  if( in == NULL ) {
118  cerr << "File not found: " << args_info.inputs[0] << endl;
119  exit(1);
120  }
121  IGBheader* head_in= new IGBheader( in );
122  if( head_in->read(true) != 0 ) gzrewind( in );
123 
124  // just output the current header
125  if( argc==2 ) {
126  output_header( head_in );
127  gzclose( in );
128  exit(0);
129  }
130 
131  IGBheader* head_out = args_info.convert_data_flag ? new IGBheader(head_in) : head_in;
132 
133  if( args_info.x_given ) head_out->x( static_cast<size_t>(args_info.x_arg) );
134  if( args_info.y_given ) head_out->y( static_cast<size_t>(args_info.y_arg) );
135  if( args_info.z_given ) head_out->z( static_cast<size_t>(args_info.z_arg) );
136  if( args_info.t_given ) head_out->t( static_cast<size_t>(args_info.t_arg) );
137  if( args_info.data_type_given ) head_out->type( args_info.data_type_arg );
138  if( args_info.system_given ) head_out->systeme( args_info.system_arg==system_arg_big ? "big_endian" : "little_endian" );
139  if( args_info.dim_x_given ) head_out->dim_x( args_info.dim_x_arg );
140  if( args_info.dim_y_given ) head_out->dim_y( args_info.dim_y_arg );
141  if( args_info.dim_z_given ) head_out->dim_z( args_info.dim_z_arg );
142  if( args_info.dim_t_given ) head_out->dim_t( args_info.dim_t_arg );
143  if( args_info.org_x_given ) head_out->org_x( args_info.org_x_arg );
144  if( args_info.org_y_given ) head_out->org_y( args_info.org_y_arg );
145  if( args_info.org_z_given ) head_out->org_z( args_info.org_z_arg );
146  if( args_info.org_t_given ) head_out->org_t( args_info.org_t_arg );
147  if( args_info.inc_x_given ) head_out->inc_x( args_info.inc_x_arg );
148  if( args_info.inc_y_given ) head_out->inc_y( args_info.inc_y_arg );
149  if( args_info.inc_z_given ) head_out->inc_z( args_info.inc_z_arg );
150  if( args_info.inc_t_given ) head_out->inc_t( args_info.inc_t_arg );
151  if( args_info.x_units_given ) head_out->unites_x( args_info.x_units_arg );
152  if( args_info.y_units_given ) head_out->unites_y( args_info.y_units_arg );
153  if( args_info.z_units_given ) head_out->unites_z( args_info.z_units_arg );
154  if( args_info.t_units_given ) head_out->unites_t( args_info.t_units_arg );
155  if( args_info.clear_comment_given ) head_out->comment(NULL);
156  if( args_info.comment_given ) head_out->comment( args_info.comment_arg );
157  if( args_info.data_factor_given ) head_out->facteur(args_info.data_factor_arg);
158  if( args_info.data_zero_given ) head_out->zero( args_info.data_zero_arg );
159  if( args_info.author_given ) head_out->aut_name( args_info.author_arg );
160  if( args_info.transparent_given ){
161  if( strlen(args_info.transparent_arg)!=2*Data_Size[head_out->type()] ){
162  cerr << "Incorrect tranparent value specified\n";
163  exit(1);
164  }
165  // convert hex digits to bytes
166  char* v = new char[Data_Size[head_out->type()]];
167  char s[3], *pp;
168  s[2] = '\0';
169  for( int i=0; i<Data_Size[head_out->type()]; i++ ){
170  s[0] = args_info.transparent_arg[i*2];
171  s[1] = args_info.transparent_arg[i*2+1];
172  v[i] = strtol( s, &pp, 16 );
173  }
174  head_out->transparent( v );
175  }
176  if( args_info.no_transparent_given ) head_out->transparent( NULL );
177 
178  // count the number of data bytes
179  if( args_info.jive_time_given )
180  jive( head_in, in );
181 
182  // make a temporary file
183  string tmpfn = ".";
184  tmpfn += basename(args_info.inputs[0]);
185  tmpfn += ".tmp";
186 
187  // determine output file
188  string ofn;
189  if( args_info.frankenstein_given ){
190  ofn = args_info.frankenstein_arg;
191  gzclose( in );
192  IGBheader *h = new IGBheader(in=gzopen( args_info.frankenstein_arg,"r" ));
193  if( h->read() != 0 ) gzrewind( in );
194  delete h;
195  }
196  else
197  ofn = args_info.inputs[0];
198 
199  if(args_info.output_file_given)
200  ofn = args_info.output_file_arg;
201 
202  FILE_PTR *out_fp;
203  if( iscompressed(ofn) ) {
204  gzFile gzout;
205  head_out->fileptr(gzout = gzopen( tmpfn.c_str(), "w" ) );
206  out_fp = new FILE_PTR( gzout );
207  } else {
208  FILE* fout;
209  head_out->fileptr(fout= fopen( tmpfn.c_str(), "w" ) );
210  out_fp = new FILE_PTR( fout );
211  }
212 
213  if( !args_info.decapitate_given )
214  head_out->write();
215 
216  if( args_info.create_given )
217  write_to_binary( head_out, in, out_fp, args_info.create_arg );
218  else if( args_info.convert_data_flag ) {
219  double *slice_data = new double[head_in->slice_sz()];
220  for( int t=0; t<head_in->t(); t++ ) {
221  head_in->read_data( slice_data );
222  head_out->write_data( slice_data );
223  }
224  delete [] slice_data;
225  } else if( args_info.transpose_flag ) {
226  int ds = head_out->data_size();
227  long offset = head_out->t()*ds;
228  long datStart = gztell( in );
229  char *buf = (char *)malloc( head_out->slice_sz()*ds );
230  for( int i=0; i<head_out->t(); i++ ) {
231  for( int j=0; j<head_out->slice_sz(); j++ ) {
232  gzseek( in,datStart+offset*j+i*ds, SEEK_SET );
233  gzread( in, buf+j*ds, ds*1 );
234  }
235  out_fp->write( buf, ds*head_out->slice_sz() );
236  }
237  delete [] buf;
238  } else {
239  const int bufsize=1000000;
240  char* buf = new char[bufsize];
241  int nb;
242  while( (nb=gzread( in, buf, bufsize )) > 0 ) out_fp->write( buf, nb );
243  delete [] buf;
244  }
245 
246  out_fp->close();
247  gzclose( in );
248  rename( tmpfn.c_str(), ofn.c_str() );
249  delete out_fp;
250  return 0;
251 }
252 
253 void output_header( IGBheader* header )
254 {
255  bool tf;
256  printf( "x dimension:\t%zu\n", header->x() );
257  printf( "y dimension:\t%zu\n", header->y() );
258  printf( "z dimension:\t%zu\n", header->z() );
259  printf( "t dimension:\t%zu\n", header->t() );
260  printf( "data type:\t%s\n", Header_Type[header->type()] );
261  header->unites(tf);
262  if( tf ) {
263  printf( "Pixel units:\t%s\n", header->unites() );
264  }
265  if( header->transparent() != NULL ) {
266  printf( "Transparent:\t%s\n", header->transparentstr() );
267  }
268  header->zero(tf);
269  if( tf ) {
270  printf( "Pixel zero:\t%g\n", header->zero() );
271  }
272  header->unites(tf);
273  if( tf ) {
274  printf( "Pixel scaling:\t%g\n", header->facteur() );
275  }
276  header->dim_x(tf);
277  if( tf ) {
278  printf( "X size:\t\t%g\n", header->dim_x() );
279  }
280  header->unites_x(tf);
281  if( tf ) {
282  printf( "X units:\t%s\n", header->unites_x() );
283  }
284  header->inc_x(tf);
285  if( tf ) {
286  printf( "Increment in x:\t%g\n", header->inc_x() );
287  }
288  header->org_x(tf);
289  if( tf ) {
290  printf( "X origin:\t%g\n", header->org_x() );
291  }
292  header->dim_y(tf);
293  if( tf ) {
294  printf( "Y size:\t\t%g\n", header->dim_y() );
295  }
296  header->unites_y(tf);
297  if( tf ) {
298  printf( "Y units:\t%s\n", header->unites_y() );
299  }
300  header->inc_y(tf);
301  if( tf ) {
302  printf( "Increment in y:\t%g\n", header->inc_y() );
303  }
304  header->org_y(tf);
305  if( tf ) {
306  printf( "Y origin:\t%g\n", header->org_y() );
307  }
308  header->dim_z(tf);
309  if( tf ) {
310  printf( "Z size:\t\t%g\n", header->dim_z() );
311  }
312  header->unites_z(tf);
313  if( tf ) {
314  printf( "Z units:\t%s\n", header->unites_z() );
315  }
316  header->inc_z(tf);
317  if( tf ) {
318  printf( "Increment in z:\t%g\n", header->inc_z() );
319  }
320  header->org_z(tf);
321  if( tf ) {
322  printf( "Z origin:\t%g\n", header->org_z() );
323  }
324  header->dim_t(tf);
325  if( tf ) {
326  printf( "T size:\t\t%g\n", header->dim_t() );
327  }
328  header->unites_t(tf);
329  if( tf ) {
330  printf( "T units:\t%s\n", header->unites_t() );
331  }
332  header->inc_t(tf);
333  if( tf ) {
334  printf( "Increment in t:\t%g\n", header->inc_t() );
335  }
336  header->org_t(tf);
337  if( tf ) {
338  printf( "T origin:\t%g\n", header->org_t() );
339  }
340  printf( "Created on:\t%s\n", header->systemestr());
341  header->aut_name(tf);
342  if( tf ) {
343  printf( "Author:\t\t%s\n", header->aut_name() );
344  }
345  if( header->comment() != NULL ){
346  char **cm = header->comment();
347  int i=0;
348  while( cm[i] != NULL )
349  printf( "#%s\n", cm[i++] );
350  }
351 }
352 
354 int jive(IGBheader *header, gzFile in )
355 {
356  z_off_t zo = gztell( in );
357  const int bufsize=8196;
358  char buff[bufsize];
359  long nb=0, nr;
360  while( (nr=gzread( in, buff, bufsize )) == bufsize )
361  nb += nr;
362  nb += nr;
363  int slicesize = header->x()*header->y()*header->z()*header->data_size();
364  header->t( nb/slicesize );
365  gzseek( in, zo, SEEK_SET );
366  if( header->dim_t() )
367  header->dim_t( (header->t()-1)*header->inc_t() );
368  return(slicesize);
369 }
370 
373 void write_slice( IGBheader *hout, IGBheader *hin, string ofnb, double *buf, int slnum )
374 {
375  long nr;
376 
377  string ofn = ofnb;
378  if( ofn.rfind(".gz")==ofn.length()-3 ) ofn.resize(ofn.length()-3);
379  if( ofn.rfind(".igb")==ofn.length()-4 ) ofn.resize(ofn.length()-3);
380 
381  char fnm[256];
382  if( (nr=hin->read_data(buf)) == hin->slice_sz() ) {
383 
384  snprintf( fnm, sizeof fnm, "%s%d.%s", ofn.c_str(), slnum, "igb" );
385 
386  if( iscompressed(ofnb) ) {
387  gzFile gzout;
388  hout->fileptr(gzout = gzopen( strcat(fnm,".gz"), "w" ) );
389  } else {
390  FILE * fout;
391  hout->fileptr(fout= fopen( fnm, "w" ) );
392  }
393 
394  hout->write();
395  hout->write_data( buf );
396  hout->close();
397  } else {
398  fprintf(stderr, "Not enough slices in file\n" );
399  exit(-1);
400  }
401 }
402 
403 bool iscompressed (string fn)
404 {
405  return fn.rfind(".gz") == fn.length()-3;
406 }
407 
int type(void)
Definition: IGBheader.h:280
void inc_x(float a)
Definition: IGBheader.h:320
void write(void *buf, int n)
Definition: igbhead.cc:50
int read_data(T *dp, size_t numt=1, char *buf=NULL)
Definition: IGBheader.h:431
bool iswrite(int i)
void dim_t(float a)
Definition: IGBheader.h:341
size_t z(void)
Definition: IGBheader.h:274
size_t slice_sz()
Definition: IGBheader.h:263
void to_bin(void *buf, T d)
Definition: IGBheader.h:529
size_t t(void)
Definition: IGBheader.h:277
void unites_x(const char *a)
Definition: IGBheader.h:353
void output_header(IGBheader *)
Definition: igbhead.cc:253
void dim_x(float a)
Definition: IGBheader.h:332
size_t data_size(void)
Definition: IGBheader.h:267
int write()
Definition: IGBheader.cc:343
void dim_y(float a)
Definition: IGBheader.h:335
void unites(const char *a)
Definition: IGBheader.h:365
void inc_y(float a)
Definition: IGBheader.h:323
int jive(IGBheader *header, gzFile in)
Definition: igbhead.cc:354
void write_slice(IGBheader *header, IGBheader *in, string ofnb, double *buf, int)
Definition: igbhead.cc:373
unsigned short num_components()
Definition: IGBheader.h:378
float org_t(void)
Definition: IGBheader.h:317
void unites_t(const char *a)
Definition: IGBheader.h:362
void close(void)
Definition: IGBheader.h:259
int read(bool quiet=false)
Definition: IGBheader.cc:757
int systeme(void)
Definition: IGBheader.h:297
void close()
Definition: igbhead.cc:49
void comment(const char *)
Definition: IGBheader.cc:1239
void unites_y(const char *a)
Definition: IGBheader.h:356
void unites_z(const char *a)
Definition: IGBheader.h:359
void aut_name(const char *a)
Definition: IGBheader.h:371
int main(int argc, char *argv[])
Definition: igbhead.cc:102
void zero(float a)
Definition: IGBheader.h:347
void facteur(float a)
Definition: IGBheader.h:344
void inc_t(float a)
Definition: IGBheader.h:329
void transparent(void *a)
Definition: IGBheader.h:374
size_t x(void)
Definition: IGBheader.h:268
void printfl(float f)
Definition: igbhead.cc:52
void inc_z(float a)
Definition: IGBheader.h:326
void write_to_binary(IGBheader *h, gzFile in, FILE_PTR *out_fp, int skip)
Definition: igbhead.cc:68
const char * Header_Type[]
Definition: IGBheader.cc:249
FILE_PTR(gzFile g)
Definition: igbhead.cc:47
void write_data(T *dp, size_t numt=1, char *buf=NULL)
Definition: IGBheader.h:397
float org_z(void)
Definition: IGBheader.h:314
size_t y(void)
Definition: IGBheader.h:271
FILE_PTR(FILE *f)
Definition: igbhead.cc:48
bool iscompressed(string fn)
Definition: igbhead.cc:403
void fileptr(gzFile f)
Definition: IGBheader.cc:329
unsigned short Data_Size[]
Definition: IGBheader.cc:264
float org_y(void)
Definition: IGBheader.h:311
char * transparentstr(void)
Definition: IGBheader.h:377
const char * systemestr(void)
Definition: IGBheader.cc:1263
float org_x(void)
Definition: IGBheader.h:308
void dim_z(float a)
Definition: IGBheader.h:338