openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
IGBextract.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 
23 #include<stdlib.h>
24 #include<stdio.h>
25 #include<iostream>
26 #include "IGBheader.h"
27 #include <assert.h>
28 #include "ext_cmdline.h"
29 #include <vector>
30 #include <algorithm>
31 
32 using namespace std;
33 
50 bool
51 process_specifier( char *sp, vector<int> &nodes, int max )
52 {
53  int inc = 1;
54  if( strchr( sp, ':' ) ){
55  char *colptr = strchr( sp, ':' );
56  int consumed;
57  sscanf( colptr+1, "%d%n", &inc, &consumed );
58  if( !consumed || consumed!=strlen(colptr+1) ) {
59  fprintf( stderr, "Illegal range specified: %s\n", sp );
60  exit(1);
61  }
62  *colptr = '\0';
63  }
64 
65  if( !strcmp( "*", sp ) ) {
66  nodes.clear();
67  for( int i=0; i<max; i+=inc )
68  nodes.push_back(i);
69  return inc==1;
70  }
71 
72  if( strchr( sp, '-' ) ){
73  int s;
74  char fs[100];
75  int nr = sscanf( sp, "%d-%s", &s, fs );
76  if( nr<2 ) {
77  fprintf( stderr, "Illegal range specified: %s\n", sp );
78  exit(1);
79  }
80 
81  int f;
82  if( !strcmp(fs, "*") )
83  f = max-1;
84  else {
85  int consumed;
86  sscanf( fs, "%d%n", &f, &consumed );
87  if( consumed != strlen(fs) ) {
88  fprintf( stderr, "Illegal range specified: %s\n", sp );
89  exit(1);
90  }
91  }
92 
93  if( s<0 || f<s ) {
94  fprintf( stderr, "Illegal range specified: %s\n", sp );
95  exit(1);
96  }
97  for( int i=s; i<=f; i+=inc )
98  nodes.push_back( i );
99 
100  } else {
101 
102  int num;
103  sscanf( sp, "%d", &num );
104 
105  if( !sscanf( sp, "%d", &num ) || num<0 ) {
106  fprintf( stderr, "Illegal node number specified in list\n" );
107  exit(1);
108  }
109  nodes.push_back( num );
110  }
111  return false;
112 }
113 
114 
115 void
116 extract_nodes( const char *list, vector<int> &nodes, int max )
117 {
118  char *lc = strdup( list );
119  char *sp = strtok( lc, "," );
120  while( sp != NULL ) {
121 
122  if( process_specifier( sp, nodes, max ) ) return;
123 
124  sp = strtok( NULL, "," );
125  }
126  free( lc );
127 }
128 
129 
130 void
131 read_nodes( char *file, vector<int> &nodes, int max )
132 {
133  FILE *in = fopen( file, "r" );
134  if( in == NULL ) {
135  fprintf( stderr, "cannot open node list file\n" );
136  exit(1);
137  }
138 #define BUFSIZE 1024
139  char buf[BUFSIZE];
140  while( fgets( buf, BUFSIZE, in ) != NULL ) {
141  if( process_specifier( buf, nodes, max ) )
142  return;
143  }
144  fclose( in );
145 }
146 
147 
156 string
157 outfile( string specified, enum_format format, int s=-1 )
158 {
159  if( specified=="-" )
160  return specified.c_str();
161 
162  string ext;
163  switch( format ) {
164  case format_arg_binary:
165  case format_arg_asciiTm:
166  case format_arg_ascii:
167  case format_arg_ascii_1cpL:
168  case format_arg_ascii_1pL:
169  case format_arg_ascii_1pLn:
170  // no standard extension
171  break;
172  case format_arg_dat_t:
173  ext = ".dat_t";
174  break;
175  case format_arg_IGB:
176  ext = ".igb";
177  break;
178  }
179 
180  string slice;
181  if( s>=0 ) {
182  char paddedS[10];
183  sprintf( paddedS, "_%05d", s );
184  slice = paddedS;
185  }
186 
187  // remove possible extension
188  if( specified.length()>ext.length() &&
189  specified.substr(specified.length()-ext.length()) == ext )
190  specified = specified.substr(0,specified.length()-ext.length());
191 
192  return specified + slice + ext;
193 }
194 
195 
207 void
208 output( vector<int> nodes, IGBheader *h, char *ofname,
209  enum_format of, int t0, int t1, int stride, bool explode, float scale)
210 {
211  IGBheader hout = *h;
212  FILE *out=NULL;
213  if( !strcmp(ofname, "-") )
214  out = stdout;
215  else if( !explode ) {
216  string out_fname = outfile(ofname, of);
217  out = fopen(out_fname.c_str(), "w");
218  if (!out) {
219  fprintf(stderr, "Can't open output file '%s'\n", out_fname.c_str());
220  exit(1);
221  }
222  }
223 
224  if (t1 >= h->t()) t1 = h->t()-1;
225  int tend = t1==-1? h->t()-1 : t1;
226 
227  hout.x( nodes.size() );
228  hout.y(1);
229  hout.z(1);
230  hout.t( explode ? 1 : (tend-t0)/stride+1 );
231  hout.org_t( h->org_t()+t0*h->inc_t() );
232  hout.inc_t(h->inc_t()*stride);
233  hout.dim_t( (t1-t0)*hout.inc_t() );
234  if( h->num_components()==3 )
235  hout.type( IGB_VEC3_f );
236  else if( h->num_components()==4 )
237  hout.type( IGB_VEC4_f );
238  else if( h->num_components()==9 )
239  hout.type( IGB_VEC9_f );
240  else
241  hout.type(IGB_FLOAT);
242 
243  if( of==format_arg_IGB ) {
244  hout.fileptr( out );
245  if( !explode ) hout.write();
246  }
247  const bool lessnodes = nodes.size() != (size_t) h->x();
248  std::vector<float> slice(h->slice_sz() * hout.num_components());
249  std::vector<float> slice_cropped(lessnodes ? nodes.size() * hout.num_components() : 0);
250  float * slice_ptr = slice.data();
251  float *output_slice = lessnodes ? slice_cropped.data() : slice_ptr;
252 
253  if( of==format_arg_dat_t )
254  fprintf( out, "%d\n", (tend-t0)/stride+1 );
255  else if( of==format_arg_ascii_1pLn )
256  fprintf( out, "%d\n", int(nodes.size()) );
257 
258  for( int t=t0; t<=tend; t+=stride ) {
259 
260  if( explode && out!=stdout ) {
261  string out_fname = outfile(ofname, of, t);
262  out = fopen(out_fname.c_str(), "w");
263  if( of==format_arg_IGB ) {
264  hout.fileptr(out);
265  hout.write();
266  }
267  }
268 
269  if( of==format_arg_dat_t )
270  fprintf(out, "%d\n", int(nodes.size()));
271 
272  h->go2slice(t);
273  h->read_data( slice_ptr );
274 
275  if( of== format_arg_asciiTm )
276  fprintf( out, "%f", h->org_t()+h->inc_t()*t );
277 
278  if(lessnodes) {
279  size_t cnt = 0;
280  for(const auto p : nodes)
281  for(size_t k=0; k < hout.num_components(); k++)
282  slice_cropped[cnt++] = slice[hout.num_components() * p + k];
283  }
284  std::for_each(output_slice, output_slice + nodes.size() * hout.num_components(), [scale](float & c){ c *= scale; });
285  switch(of) {
286  case format_arg_IGB:
287  case format_arg_binary:
288  fwrite(output_slice, sizeof(float), hout.slice_sz() * hout.num_components(), out );
289  break;
290  case format_arg_asciiTm:
291  {
292  size_t cnt=0;
293  for(size_t i = 0; i < nodes.size(); i++)
294  for(size_t k=0; k < hout.num_components(); k++)
295  fprintf(out, " %g", output_slice[cnt++]);
296  break;
297  }
298  case format_arg_ascii_1cpL:
299  {
300  size_t cnt=0;
301  for(size_t i = 0; i < nodes.size(); i++)
302  for(size_t k=0; k < hout.num_components(); k++)
303  fprintf(out, " %g\n", output_slice[cnt++]);
304  break;
305  }
306  case format_arg_ascii_1pL:
307  case format_arg_ascii_1pLn:
308  {
309  size_t cnt=0;
310  for(size_t i = 0; i < nodes.size(); i++) {
311  for(size_t k=0; k < hout.num_components(); k++)
312  fprintf(out, " %g", output_slice[cnt++]);
313  fputs("\n", out);
314  }
315  break;
316  }
317  case format_arg_ascii:
318  {
319  size_t cnt=0;
320  for(size_t i=0; i < nodes.size(); i++)
321  for(size_t k=0; k < hout.num_components(); k++)
322  fprintf(out, "%g ", output_slice[cnt++]);
323  break;
324  }
325  case format_arg_dat_t:
326  {
327  size_t cnt=0;
328  for(size_t i=0; i < nodes.size(); i++)
329  for(size_t k=0; k < hout.num_components(); k++)
330  fprintf(out, "%g\n", output_slice[cnt++]);
331  break;
332  }
333  }
334  if( of==format_arg_ascii || of==format_arg_asciiTm )
335  fprintf( out, "\n" );
336  if( explode && out!=stdout )
337  fclose(out);
338  }
339 
340  if( !explode && out != stdout )
341  fclose( out );
342 }
343 
344 
345 int
346 main( int argc, char *argv[] )
347 {
348  gengetopt_args_info args;
349 
350  // let's call our cmdline parser
351  if (cmdline_parser (argc, argv, &args) != 0)
352  exit(1);
353 
354  // the last arg must be a file name
355  gzFile in = ( !args.inputs_num || !strcmp(args.inputs[0],"-") ) ?
356  gzdopen( dup(0), "r" ): // duplicate stdin
357  gzopen( args.inputs[0], "r" );
358  if( in == NULL ) {
359  cerr << "File not found: " << args.inputs[0] << endl;
360  exit(1);
361  }
362  IGBheader header(in);
363  header.read(true);
364 
365  if( args.explode_flag && args.format_arg==format_arg_dat_t ) {
366  cerr << "Cannot explode into dat_t format" << endl;
367  exit(1);
368  }
369 
370  vector<int> nodes;
371  if( args.list_given )
372  extract_nodes( args.list_arg, nodes, header.slice_sz() );
373  else if( args.input_file_given )
374  read_nodes( args.input_file_arg, nodes, header.slice_sz() );
375  else // read them all
376  extract_nodes( "*", nodes, header.slice_sz() );
377 
378  for( auto p : nodes )
379  if( header.slice_sz() <= p ) {
380  cerr << "Error: out of range nodes requested" << endl;
381  exit(1);
382  }
383 
384  int f0=args.f0_arg, f1=args.f1_arg;
385  if( args.t0_given )
386  f0 = std::round((args.t0_arg-header.org_t())/header.inc_t());
387  if( f0<0 || f0>=header.t() ) {
388  cerr << "Illegal ending frame specified: "<<f0<<" not in [0,"<<header.t()<<")"<< endl;
389  exit(1);
390  }
391  if( args.t1_given )
392  f1 = args.t1_arg==-1 ? header.t()-1 : std::round((args.t1_arg-header.org_t())/header.inc_t());
393  if( f1>=0 && f1<f0 ) {
394  cerr << "Illegal ending frame specified: "<<f1<<" not in ["<<f0<<","<<header.t()<<")"<< endl;
395  exit(1);
396  }
397 
398  output(nodes, &header, args.output_file_arg,args.format_arg,
399  f0, f1, args.df_arg, args.explode_flag, args.scale_arg );
400 
401  cmdline_parser_free(&args);
402 }
int type(void)
Definition: IGBheader.h:280
void read_nodes(char *file, vector< int > &nodes, int max)
Definition: IGBextract.cc:131
int read_data(T *dp, size_t numt=1, char *buf=NULL)
Definition: IGBheader.h:431
void dim_t(float a)
Definition: IGBheader.h:341
#define IGB_FLOAT
Definition: IGBheader.h:60
constexpr T max(T a, T b)
Definition: ion_type.h:31
size_t z(void)
Definition: IGBheader.h:274
#define IGB_VEC9_f
Definition: IGBheader.h:76
void extract_nodes(const char *list, vector< int > &nodes, int max)
Definition: IGBextract.cc:116
size_t slice_sz()
Definition: IGBheader.h:263
int go2slice(size_t s)
go to a slice
Definition: IGBheader.cc:1459
size_t t(void)
Definition: IGBheader.h:277
#define IGB_VEC3_f
Definition: IGBheader.h:71
int write()
Definition: IGBheader.cc:343
#define BUFSIZE
void output(vector< int > nodes, IGBheader *h, char *ofname, enum_format of, int t0, int t1, int stride, bool explode, float scale)
Definition: IGBextract.cc:208
unsigned short num_components()
Definition: IGBheader.h:378
float org_t(void)
Definition: IGBheader.h:317
int read(bool quiet=false)
Definition: IGBheader.cc:757
bool process_specifier(char *sp, vector< int > &nodes, int max)
Definition: IGBextract.cc:51
string outfile(string specified, enum_format format, int s=-1)
Definition: IGBextract.cc:157
void inc_t(float a)
Definition: IGBheader.h:329
size_t x(void)
Definition: IGBheader.h:268
int main(int argc, char *argv[])
Definition: IGBextract.cc:346
size_t y(void)
Definition: IGBheader.h:271
void fileptr(gzFile f)
Definition: IGBheader.cc:329
#define IGB_VEC4_f
Definition: IGBheader.h:73