openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
IGBops.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 
53 #include<stdlib.h>
54 #include<stdio.h>
55 #include<unistd.h>
56 #include<iostream>
57 #include "IGBheader.h"
58 #include <assert.h>
59 #include "ops_cmdline.h"
60 #include <vector>
61 #include "muParser.h"
62 #ifdef _OPENMP
63 #include <omp.h>
64 #endif
65 
66 const int MAX_NC=4;
67 
68 using namespace std;
69 
70 enum Data_file_t { NONE=0, SPACE_TIME=5, SPACE=4, PT_TIME=3, PT=2, TIME=1 };
71 int Tbit=1; // time bit
72 
76 class Interp {
77  public:
78  Interp( char *ind, char *weight );
79  int size() { return ind.size(); }
80  vector<float*>w;
81  vector<int*>ind;
82  int n=0;
83 };
84 
91 Interp::Interp( char *indf, char *wf )
92 {
93  char buf[8192];
94  FILE *iin = fopen( indf, "r" );
95  while( fgets( buf, 8192, iin )!=NULL ) {
96  if( !n ) {
97  char *l1 = strdup(buf);
98  char *ptr = strtok( l1, " " );
99  while((ptr=strtok( NULL, " " ))) n++;
100  free(l1);
101  }
102  int *ip = new int[n];
103  istringstream oss( buf );
104  for( int i=0; i<n; i++ )
105  oss>>ip[i];
106  ind.push_back(ip);
107  }
108  FILE *win = fopen( wf, "r" );
109  for( int i=0; i<ind.size(); i++ ) {
110  fgets( buf, 8192, win );
111  float *wp = new float[n];
112  istringstream oss( buf );
113  for( int j=0; j<n; j++ )
114  oss >> wp[j];
115  w.push_back( wp );
116  }
117 }
118 
119 
120 
121 #define data_t float
123 
124 void output_scratch( data_t* s, FILE *out, int n, bool bin );
125 void output_scratch( vector<data_t>& s, FILE *out, bool bin );
126 
127 
128 bool overlap(IGBheader *first, IGBheader *second )
129 {
130  return fabs(second->org_t() - (first->org_t()+(first->t()-1)*first->inc_t()) ) < first->inc_t()/2.;
131 }
132 
133 
134 // concatenate X and Y in time
135 void joinXY( IGBheader *h1, IGBheader *h2, IGBheader *hout )
136 {
137  if( h2==NULL ) {
138  fprintf( stderr, "Y must be an IGB file\n\n" );
139  exit(1);
140  }
141 
142  gzFile in1 = (gzFile)(h1->fileptr());
143  gzFile in2 = (gzFile)(h2->fileptr());
144  FILE* out = (FILE *)(hout->fileptr());
145 
146  std::vector<data_t> dat1(hout->slice_sz() * hout->num_components());
147  for( int i=0; i<h1->t(); i++ ) {
148  h1->read_data( dat1.data() );
149  fwrite( dat1.data(), sizeof(data_t), hout->slice_sz() * hout->num_components(), out );
150  }
151 
152  bool skip_second_torg = overlap( h1, h2 );
153  for( int i=0; i<h2->t(); i++ ) {
154  h2->read_data( dat1.data() );
155  if(skip_second_torg)
156  skip_second_torg = false;
157  else
158  fwrite( dat1.data(), sizeof(data_t), hout->slice_sz() * hout->num_components(), out );
159  }
160  fclose(out);
161 }
162 
163 // put Y beside X
164 void join_space( data_t *x, int nx, data_t *y, int ny, FILE *out )
165 {
166  fwrite( x, sizeof(data_t), nx, out );
167  fwrite( y, sizeof(data_t), ny, out );
168 }
169 
170 
171 // evaluate arbitrary expression for one file
172 void eval_expr( data_t *x, int n, char *expr, FILE *out, data_t *result )
173 {
174  using namespace mu;
175 
176  static Parser hParser;
177  static bool init=false;
178 
179  if( !init ) {
180  hParser.DefineVar( "X", x );
181  hParser.SetExpr(expr);
182  init = true;
183  }
184  try {
185  hParser.Eval(result, n);
186  }
187  catch (Parser::exception_type &e)
188  {
189  std::cerr << e.GetMsg() << std::endl;
190  exit(1);
191  }
192  output_scratch( result, out, n, true );
193 }
194 
195 
196 // evaluate arbitrary expression as function of two files
197 // nc>0 means that x1 is a point
208 template<class T>
209 void eval_expr( data_t *x0, T* x1, int n, char *expr, FILE *out, int nc=0 )
210 {
211  using namespace mu;
212 
213  Parser hParser;
214  data_t a, b;
215  hParser.DefineVar( "X", &a );
216  hParser.DefineVar( "Y", &b );
217  hParser.SetExpr(expr);
218 
219  data_t *y = new data_t[n];
220 
221  for( int i=0; i<n; i++ ) {
222  try {
223  a = x0[i];
224  b = x1[nc>0 ? i%nc : i];
225  y[i] = hParser.Eval();
226  }
227  catch (Parser::exception_type &e)
228  {
229  std::cerr << e.GetMsg() << std::endl;
230  exit(1);
231  }
232  }
233  output_scratch( y, out, n, true );
234  delete[] y;
235 }
236 
237 
238 // central difference
239 void center_diff( data_t* x, int slsz, vector<data_t>& odat, int t, float dt, FILE *out)
240 {
241  int older= (t%2) ? 0 : slsz;
242  int old = older ? 0 : slsz;
243 
244  if( t==1 ) { // forward difference for first frame
245  for(int i=0;i<slsz;i++){
246  data_t d = (x[i]-odat[old+i])/dt;
247  fwrite( &d, sizeof(data_t), 1, out );
248  }
249  } else if( t>1 ) {
250  for(int i=0;i<slsz;i++){
251  data_t d = (x[i]-odat[older+i])/dt/2.;
252  fwrite( &d, sizeof(data_t), 1, out );
253  }
254  }
255  for(int i=0;i<slsz;i++)
256  odat[older+i] = x[i];
257 }
258 
259 
260 // forward difference
261 void forward_diff( data_t* x, int slsz, vector<data_t>& old, int t, float dt, FILE *out)
262 {
263  for(int i=0;i<slsz;i++){
264  if( t ) {
265  data_t d = (x[i]-old[i])/dt;
266  fwrite( &d, sizeof(data_t), 1, out );
267  }
268  old[i] = x[i];
269  }
270 }
271 
272 
273 // aX[i] + bY[i%nc]
274 void aXPLUS_bY( data_t* x, data_t *y, data_t a, data_t b, int slsz, int nc, FILE *out)
275 {
276  for(int i=0;i<slsz;i++){
277  data_t datum = x[i]*a + b*y[i%nc];
278  fwrite( &datum, sizeof(data_t), 1, out );
279  }
280 }
281 
282 // aX[i] + bY[i]
283 template<class T>
284 void aXPLUS_bY( data_t* x, T* y, data_t a, data_t b, int slsz, FILE *out)
285 {
286  for(int i=0;i<slsz;i++){
287  data_t datum = x[i]*a + b*y[i];
288  fwrite( &datum, sizeof(data_t), 1, out );
289  }
290 }
291 
292 // avgX[*,t]
293 void avgXt( data_t* x, int slsz, int nc, float t, FILE *out)
294 {
295  data_t sumf[MAX_NC];
296  for( int c=0; c<nc; c++ ) {
297  double s=0;
298 #pragma omp parallel for reduction( +:s )
299  for(int i=c;i<slsz;i+=nc)
300  s += x[i];
301  sumf[c] = s/double(slsz);
302  }
303  fwrite( &sumf, sizeof(sumf[0]), nc, out );
304 }
305 
306 
307 // avgX[i,*]
308 void avgXi( data_t* x, vector<data_t>& sum, const data_t fac)
309 {
310 #pragma omp parallel for
311  for(int i=0; i<sum.size(); i++){
312  sum[i] += fac * x[i];
313  }
314 }
315 
316 void varXi( data_t* x, data_t* mean, vector<data_t>& sum, const data_t fac)
317 {
318 #pragma omp parallel for
319  for(int i=0; i<sum.size(); i++){
320  sum[i] += fac * (x[i] - mean[i]) * (x[i] - mean[i]);
321  }
322 }
323 
324 // aX[i,t]+b
325 void aXPLUS_b( data_t* x, data_t a, data_t b, int slsz, FILE *out)
326 {
327  for(int i=0;i<slsz;i++){
328  data_t datum = x[i]*a+b;
329  fwrite( &datum, sizeof(data_t), 1, out );
330  }
331 }
332 
333 // X[i]/Y[i]
334 void XdivY( data_t* x, data_t* y, int slsz, FILE *out)
335 {
336  bool NAN_FLAG = 0;
337  bool INF_FLAG = 0;
338 
339  for(int i=0;i<slsz;i++){
340  data_t datum = x[i]/y[i];
341  if (std::isinf(datum)) {
342  datum = 0.; // handle infinite numbers
343  INF_FLAG = 1;
344  }
345  if (std::isnan(datum)){
346  datum = 0.; // handle not a number
347  NAN_FLAG = 1;
348  }
349  fwrite( &datum, sizeof(data_t), 1, out );
350  }
351  if (INF_FLAG) std::cout << "Infinite number(s) detected (and set to zero)" << endl;
352  if (NAN_FLAG) std::cout << "NaN(s) detected (and set to zero)" << endl;
353 }
354 
355 
356 // X[i]/Y[0]
357 void XdivY( data_t* x, data_t* y, int slsz, int nc, FILE *out)
358 {
359  bool NAN_FLAG = 0;
360  bool INF_FLAG = 0;
361 
362  for(int i=0;i<slsz;i++){
363  data_t datum = x[i]/y[i%nc];
364  if (std::isinf(datum)) {
365  datum = 0.; // handle infinite numbers
366  INF_FLAG = 1;
367  }
368  if (std::isnan(datum)){
369  datum = 0.; // handle not a number
370  NAN_FLAG = 1;
371  }
372  fwrite( &datum, sizeof(datum), 1, out );
373  }
374  if (INF_FLAG) std::cout << "Infinite number(s) detected (and set to zero)" << endl;
375  if (NAN_FLAG) std::cout << "NaN(s) detected (and set to zero)" << endl;
376 }
377 
378 // X[i]*Y[i]
379 void XmulY( data_t* x, data_t* y, int slsz, FILE *out)
380 {
381  for(int i=0;i<slsz;i++){
382  data_t datum = x[i]*y[i];
383  fwrite( &datum, sizeof(data_t), 1, out );
384  }
385 }
386 
387 // X[i]*Y[0]
388 void XmulY( data_t* x, data_t* y, int slsz, int nc, FILE *out)
389 {
390  for(int i=0;i<slsz;i++){
391  data_t datum = x[i]*y[i%nc];
392  fwrite( &datum, sizeof(datum), 1, out );
393  }
394 }
395 
396 // MIN_i( X[i,t] )
397 void minXt( data_t* x, int slsz, int nc, float t, FILE *out )
398 {
399  data_t min[MAX_NC];
400  for( int i=0; i<nc; i++ ) min[i]=x[i];
401 #pragma omp parallel for
402  for(int i=nc;i<slsz*nc;i++){
403  if( x[i]<min[i%nc] )
404  min[i%nc] = x[i];
405  }
406  fwrite( min, sizeof(data_t), nc, out );
407 }
408 
409 
410 // t(MAX_t( X[i,t] ))
411 void t_maxX( data_t* x, data_t* mx, data_t* tm, size_t n, float t )
412 {
413 #pragma omp parallel for
414  for( int i=0; i<n; i++ ){
415  if( x[i]>mx[i] ){
416  mx[i] = x[i];
417  tm[i] = t;
418  }
419  }
420 }
421 
422 // MAX_t( X[i,t] )
423 void maxXi( data_t* x, data_t* mx, size_t n )
424 {
425 #pragma omp parallel for
426  for( int i=0; i<n; i++ ){
427  if( x[i]>mx[i] )
428  mx[i] = x[i];
429  }
430 }
431 
432 
433 // MIN_t( X[i,t] )
434 void minXi( data_t* x, data_t* mn, size_t n )
435 {
436 #pragma omp parallel for
437  for( int i=0; i<n; i++ ){
438  if( x[i]<mn[i] )
439  mn[i] = x[i];
440  }
441 }
442 
443 
444 // MAX_i( X[i,t] )
445 void maxXt( data_t* x, int slsz, int nc, float t, FILE *out )
446 {
447  data_t max[MAX_NC];
448  for( int i=0; i<nc; i++ ) max[i]=x[i];
449 #pragma omp parallel for
450  for(int i=nc;i<slsz*nc;i++){
451  if( x[i]>max[i%nc] )
452  max[i%nc] = x[i];
453  }
454  fwrite( max, sizeof(data_t), nc, out );
455 }
456 
457 // X[i,t]+a
458 void XPLUS_a( data_t* x, data_t a, int slsz, FILE *out )
459 {
460  for(int i=0;i<slsz;i++){
461  data_t datum = x[i]+a;
462  fwrite( &datum, sizeof(data_t), 1, out );
463  }
464 }
465 
466 
467 // apply a simple box filter to the xXy grid where the center pixel
468 // is the average of the enclosing dXd square, do not filter edges
469 void Xbox( data_t *dat, int x, int y, FILE *out, int d, int nc )
470 {
471  assert(nc==1);
472 
473  int fm = d/2,
474  fp = d - 1 - fm;
475  float area = d*d;
476  for( int j=0; j<y; j++ )
477  for( int i=0; i<x; i++ ) {
478  data_t datum = 0;
479  if( j<fm || j>y-fp-1 || i<fm || i>x-fp-1 )
480  datum = dat[j*x+i];
481  else {
482  for( int m=j-fm; m<=j+fp; m++ )
483  for( int n=i-fm; n<=i+fp; n++ )
484  datum += dat[m*x+n];
485  datum /= area;
486  }
487  fwrite( &datum, sizeof(data_t), 1, out );
488  }
489 }
490 
491 
492 // soft clipping: remove offending point and interpolate
493 void sclip_ab( data_t *dat, int sz, data_t upper, data_t lower, vector<data_t>& scr, int t,
494  int tfinal, FILE *out )
495 {
496  unsigned int curr = (t%3)*sz;
497  unsigned int old = (curr+2*sz)%(3*sz);
498  unsigned int older = (curr+sz)%(3*sz);
499 
500  memcpy( scr.data()+curr, dat, sizeof(data_t)*sz );
501 
502  if( !t )
503  fwrite( dat, sizeof(data_t), sz, out );
504 
505  if( t>1 ) {
506  for( int i=0; i<sz; i++ )
507  if( scr[old+i]>upper || scr[old+i]<lower ) {
508  if( scr[curr+i]>upper || scr[curr+i]<lower )
509  scr[old+i] = scr[older+i];
510  else
511  scr[old+i] = (scr[older+i]+scr[curr+i])/2.;
512  }
513  fwrite( scr.data()+old, sizeof(data_t), sz, out );
514  }
515 
516  if( t==tfinal-1 )
517  fwrite( dat, sizeof(data_t), sz, out );
518 }
519 
520 
528 void
529 magnitude( data_t* x, int n, int nc, FILE *out)
530 {
531  for(int i=0; i<n; i+=nc){
532  data_t mag = 0;
533  for( int j=0; j<nc; j++ )
534  mag += x[i+j]*x[i+j];
535  mag = sqrt( mag );
536  fwrite( &mag, sizeof(mag), 1, out );
537  }
538 }
539 
549 void
550 dot( data_t* x, data_t *y, int n, int nc, FILE *out, bool pt)
551 {
552  int k = !pt;
553  for(int i=0; i<n; i+=nc){
554  data_t dot = 0;
555  for( int j=0; j<nc; j++ )
556  dot += x[i+j]*y[k*i+j];
557  fwrite( &dot, sizeof(dot), 1, out );
558  }
559 }
560 
570 template<class T>
571 void
572 simp38( T *x, int n, int t, data_t *runInt, float dt )
573 {
574  int s=t%3;
575  T scale = dt;
576 
577  if( t<0 )
578  scale *= 3./8.;
579  else if( s==0 )
580  scale *= 3./4.;
581  else if( s==1 || s==2 )
582  scale *= 9./8.;
583 
584 #pragma omp parallel for
585  for( int i=0; i<n; i++ )
586  runInt[i] += x[i]*scale;
587 }
588 
589 
592 template<class T>
593 void
594 midpt( T *x, int n, bool end, data_t *reduce, float dt )
595 {
596 #pragma omp parallel for
597  for( int i=0; i<n; i++ )
598  if( end )
599  reduce[i] += dt*x[i]/2.;
600  else
601  reduce[i] += dt*x[i];
602 }
603 
604 
615 void
616 interpolate( data_t *dat, vector<int*>&idx, vector<float*>& w, int n, FILE *out )
617 {
618  for( int i=0; i<idx.size(); i++ ) {
619  data_t datum=0.;
620  for( int j=0; j<n; j++ )
621  datum += dat[idx[i][j]]*w[i][j];
622  fwrite( &datum, sizeof(data_t), 1, out );
623  }
624 }
625 
626 
635 template<class T>
636 void
637 normalize( IGBheader *h, FILE *out, T *min, T* max )
638 {
639  size_t slice_mem = h->slice_sz()*h->num_components();
640  h->go2slice(0);
641  vector<data_t> slice( slice_mem );
642 
643  for( int t=0; t< h->t(); t++ ) {
644  h->read_data( slice.data() );
645  for( int j=0; j<slice_mem; j++ ){
646  data_t datum = (slice[j]-min[j])/(max[j]-min[j]);
647  fwrite( &datum, sizeof(data_t), 1, out );
648  }
649  }
650 }
651 
652 
667 void
668 maxP2P( data_t* dat, data_t* big, data_t* maxpp, int sz, int t )
669 {
670  data_t *old = big;
671  data_t *change = big+sz;
672  data_t *extrem = big+2*sz;
673 
674  if( !t ) memcpy( old, dat, sizeof(data_t)*sz );
675 
676 #pragma omp parallel for
677  for( int i=0; i<sz; i++ ) {
678  data_t dV = dat[i] - old[i];
679  if( dV*change[i] < 0 ) {
680  data_t pp = fabs( extrem[i]-old[i] );
681  if( pp > maxpp[i] ) maxpp[i] = pp;
682  change[i] = copysign( 1., dV );
683  extrem[i] = old[i];
684  } else if( change[i]==0. && dV != 0 ) {
685  change[i] = copysign( 1., dV );
686  }
687  }
688  memcpy( old, dat, sizeof(data_t)*sz );
689 }
690 
691 
692 //output scratch vector
693 void
694 output_scratch( data_t* s, FILE *out, int n, bool bin )
695 {
696  for( int i=0; i<n; i++ )
697  if( bin ) {
698  fwrite( s+i, sizeof(data_t), 1, out );
699  }else
700  fprintf( out, "%g\n", s[i] );
701 }
702 
703 //output scratch vector
704 void
705 output_scratch( vector<data_t>& s, FILE *out, bool bin=false )
706 {
707  for( auto i : s )
708  if( bin ) {
709  data_t datum = i;
710  fwrite( &datum, sizeof(data_t), 1, out );
711  }else
712  fprintf( out, "%g\n", i );
713 }
714 
715 
716 
717 // temporal running average
718 void ravgX_ta( data_t *dat1, int sz, int t, int a, vector<data_t>& tmp, FILE *out )
719 {
720  vector <data_t> sum( sz, 0 );
721 
722 #pragma omp parallel for
723  for( int i=0; i<sz; i++ ) {
724 
725  tmp[i+(t%a)*sz] = dat1[i];
726 
727  for( int j=0; j<a; j++ ) {
728  sum[i] += tmp[i+j*sz];
729  }
730  sum[i] /= (float)a;
731  }
732 
733  if( t<a-1 )
734  return;
735 
736  output_scratch( sum, out, true );
737 }
738 
739 
740 //read in an ascii vector file
741 void
742 read_vector( vector<data_t> &y, gzFile in, int nc )
743 {
744  assert( nc<=MAX_NC );
745  char buf[1024];
746  float datum[MAX_NC];
747  while( gzgets( in, buf, 1024 ) != Z_NULL &&
748  sscanf( buf, "%f %f %f %f", datum, datum+1, datum+2, datum+3 )==nc )
749  for( int i=0; i<nc; i++ ) y.push_back( datum[i] );
750 }
751 
752 
758 std::string
759 cmdline_str( int n, char *w[] )
760 {
761  std::string cl;
762  for( int i=0; i<n; i++ ) {
763  cl += w[i];
764  cl += " ";
765  }
766  return cl;
767 }
768 
769 
770 int
771 main( int argc, char *argv[] )
772 {
773 #ifdef _OPENMP
774  if( !getenv("OMP_NUM_THREADS") )
775  omp_set_num_threads(1);
776 #endif
777 
778  gengetopt_args_info opts;
779 
780  std::string cmdline = cmdline_str( argc, argv );
781 
782  // let's call our cmdline parser
783  if (cmdline_parser (argc, argv, &opts) != 0)
784  exit(1);
785 
786  // the last arg must be a file name
787  gzFile in1 = (!opts.inputs_num || !strcmp(opts.inputs[0],"-")) ?
788  gzdopen( dup(0), "r" ): // duplicate stdin
789  gzopen( opts.inputs[0], "r" );
790  if( in1 == NULL ) {
791  cerr << "File not found: " << opts.inputs[0] << endl;
792  exit(1);
793  }
794  // make sure we can open X which must be an IGB file
795  IGBheader* h1= new IGBheader( in1 );
796  if( h1->read(true) & 0xff ) {
797  cerr << "File 1 is not a proper IGB file: "<< opts.inputs[0] << endl;
798  exit(1);
799  }
800 
801  Interp *interp=NULL;
802  if( opts.interpolate_flag )
803  interp = new Interp( opts.inputs[1], opts.inputs[2] );
804 
805  // is a second file necessary?
806  int numf=1;
807  if( opts.op_given )
808  switch( opts.op_arg ) {
809  case op_arg_aXPLUS_bY:
810  case op_arg_XdivY:
811  case op_arg_XmulY:
812  case op_arg_joinXY_t:
813  case op_arg_joinXY_i:
814  case op_arg_dotXY:
815  numf = 2;
816  break;
817  }
818  else if( opts.expression_given )
819  if( strchr( opts.expression_arg, 'Y' ) )
820  numf =2;
821 
822  // read header and determine size of second file if necessary
823  gzFile in2;
824  IGBheader *h2 = NULL;
825  vector<float> Y;
826  Data_file_t Ytype = NONE;
827  bool Yigb=false;
828 
829  if( numf==2 ) {
830  if( opts.inputs_num<2 ) {
831  cerr << "A 2nd file is needed : " << endl;
832  exit(1);
833  }
834  if( (in2=gzopen( opts.inputs[1], "r" ))==NULL ) {
835  cerr << "file 2 not found : " << opts.inputs[1] << endl;
836  exit(1);
837  }
838  // check file type
839  h2 = new IGBheader( in2 );
840  if( h2->read(true) & 0xff ) {
841  gzrewind( in2 );
842  read_vector( Y, in2, h1->num_components() );
843  if( Y.size()==h1->slice_sz()*h1->num_components() )
844  Ytype = SPACE;
845  else if( Y.size()==h1->t()*h1->num_components() )
846  Ytype = PT_TIME;
847  else if( Y.size()==h1->num_components() )
848  Ytype = PT;
849  else if( opts.op_arg != op_arg_joinXY_i ) {
850  cerr << "incompatile file sizes! Y must have " << h1->t()*h1->num_components()
851  << " or " << h1->slice_sz()*h1->num_components() << " or " << h1->num_components()
852  << " entries" << endl;
853  exit(1);
854  }
855  } else {
856  Yigb = true;
857  if( h2->num_components() != h1->num_components() ) {
858  cerr << "Number of data components do not match" << endl;
859  exit(1);
860  }
861  if( h2->slice_sz()==h1->slice_sz() && h2->t()==h1->t() )
862  Ytype = SPACE_TIME;
863  else if( h2->slice_sz()==h1->slice_sz() ) {
864  Ytype = SPACE;
865  if( h2->t()!=1 && opts.op_arg != op_arg_joinXY_t ) {
866  cerr << "Second IGB file time does not match" << endl;
867  exit(1);
868  }
869  Y.resize(h1->slice_sz()*h1->num_components());
870  } else if( h2->slice_sz()==1 && h2->t()==h1->t() )
871  Ytype = PT_TIME;
872  else if ( h2->t()==h1->t() )
873  Ytype = TIME;
874  else if ( h2->slice_sz()==1 && h2->t()==1 )
875  Ytype = PT;
876  else if( opts.op_arg != op_arg_joinXY_i ) {
877  cerr << "incompatile IGB file sizes!" << endl;
878  exit(1);
879  }
880  }
881  }
882 
883  // determine the output file size
884  IGBheader *hout;
885  FILE *out;
886  string fname = opts.output_file_arg;
887  if( fname=="-" )
888  out= stdout;
889  else{
890  if( fname.length()<5 || fname.substr(fname.length()-4,4)!=".igb" )
891  fname += opts.extension_arg;
892  out = fopen( fname.c_str(), "w" );
893  }
894  hout = new IGBheader;
895  *hout = *h1;
896  hout->type(IGB_DATA_T[h1->num_components()]);
897  switch( opts.op_arg ) {
898  case op_arg_ravgX_ta:
899  hout->t(h1->t()-opts.a_arg+1);
900  break;
901  case op_arg_joinXY_i:
902  if( Yigb )
903  hout->x(h1->slice_sz()+h2->slice_sz());
904  else
905  hout->x(h1->slice_sz()+Y.size());
906  hout->y(1);
907  hout->z(1);
908  break;
909  case op_arg_joinXY_t:
910  hout->type(h1->type());
911  hout->t(h1->t()+h2->t()-overlap(h1,h2));
912  break;
913  case op_arg_diff_f:
914  case op_arg_diff_c:
915  hout->t(h1->t()-1);
916  break;
917  case op_arg_minX_i:
918  case op_arg_maxX_i:
919  case op_arg_t_maxX:
920  case op_arg_avgX_i:
921  case op_arg_stdD_i:
922  case op_arg_intX_i:
923  case op_arg_maxP2P:
924  hout->t(1);
925  hout->type(IGB_DATA_T[1]);
926  break;
927  case op_arg_minX_t:
928  case op_arg_maxX_t:
929  case op_arg_avgX_t:
930  case op_arg_stdD_t:
931  hout->x(1);
932  hout->y(1);
933  hout->z(1);
934  break;
935  case op_arg_magX:
936  case op_arg_dotXY:
937  hout->type(IGB_DATA_T[1]);
938  break;
939  }
940  if( opts.interpolate_flag ) {
941  hout->x( interp->size() );
942  hout->y(1);
943  hout->z(1);
944  }
945  //hout->comment( cmdline.c_str() );
946  hout->fileptr( out );
947  hout->write();
948 
949  if( opts.op_arg == op_arg_joinXY_t ) {
950  joinXY( h1, h2, hout );
951  exit(0);
952  }
953 
954  // allocate some input memory
955  std::vector<float> dat1(h1->slice_sz() * h1->num_components());
956  std::vector<std::vector<float>> dat3;
957  if( Yigb )
958  Y.assign(h2->slice_sz()*h2->num_components(), 0.0f);
959 
960  // initialize if we need scratch data
961  int h1_slice_comps = h1->slice_sz()*h1->num_components();
962  std::vector<data_t> t_reduce; // automatically output at the end for a temporal reduction
963  std::vector<data_t> big_scratch;
964  switch( opts.op_arg ) {
965  case op_arg_minX_i:
966  case op_arg_maxX_i:
967  case op_arg_avgX_i:
968  case op_arg_stdD_i:
969  case op_arg_intX_i :
970  case op_arg_intX_t :
971  t_reduce.assign( h1_slice_comps, 0 );
972  break;
973  case op_arg_t_maxX:
974  t_reduce.assign( h1_slice_comps, 0 );
975  big_scratch.assign( h1_slice_comps, 0 );
976  break;
977  case op_arg_ravgX_ta :
978  big_scratch.assign( h1_slice_comps*opts.a_arg, 0 );
979  break;
980  case op_arg_sclip_ab :
981  big_scratch.assign( h1_slice_comps*3, 0 );
982  break;
983  case op_arg_diff_f :
984  big_scratch.assign( h1_slice_comps, 0 );
985  break;
986  case op_arg_diff_c :
987  case op_arg_norm_i :
988  big_scratch.assign( h1_slice_comps*3, 0 );
989  break;
990  case op_arg_maxP2P :
991  big_scratch.assign( h1_slice_comps*3, 0 );
992  t_reduce.assign( h1_slice_comps, 0 );
993  break;
994  }
995  if( opts.expression_given )
996  big_scratch.assign( h1_slice_comps*opts.a_arg, 0 );
997 
998  // loop through the data
999  for( int t=0; t<h1->t(); t++ ) {
1000 
1001  float tm = h1->org_t() + t*h1->inc_t();
1002 
1003  h1->read_data( dat1.data() );
1004 
1005  if(opts.op_arg == op_arg_stdD_i)
1006  dat3.push_back(dat1);
1007 
1008  if( Yigb && (Ytype&Tbit || !t) )
1009  h2->read_data( Y.data() );
1010 
1011  float *Yptr = Y.data();
1012  if( !Yigb && Ytype==PT_TIME )
1013  Yptr += h1->num_components()*t;
1014 
1015  // first data read initialization
1016  if( !t ) {
1017  if( opts.op_arg==op_arg_minX_i || opts.op_arg==op_arg_maxX_i ) {
1018  t_reduce = dat1;
1019  } else if( opts.op_arg==op_arg_norm_i ) {
1020  copy( dat1.begin(), dat1.end(), big_scratch.begin() );
1021  copy( dat1.begin(), dat1.end(), big_scratch.begin()+h1_slice_comps );
1022  } else if( opts.op_arg==op_arg_t_maxX) {
1023  big_scratch = dat1;
1024  }
1025  }
1026 
1027  if( opts.expression_given ){
1028  if( Ytype==NONE )
1029  eval_expr( dat1.data(), h1_slice_comps, opts.expression_arg, out, big_scratch.data() );
1030  else {
1031  bool space = Ytype==SPACE_TIME || Ytype==SPACE;
1032  eval_expr( dat1.data(), Yptr, h1_slice_comps, opts.expression_arg,
1033  out, space? 0 : h1->num_components() );
1034  }
1035  }
1036 
1037  if( opts.interpolate_given )
1038  interpolate( dat1.data(), interp->ind, interp->w, interp->n, out );
1039 
1040  switch( opts.op_arg ) {
1041  case op_arg_aXPLUS_bY:
1042  if( Ytype==SPACE_TIME || Ytype==SPACE )
1043  aXPLUS_bY(dat1.data(), Yptr, opts.a_arg, opts.b_arg,
1044  h1_slice_comps,out);
1045  else
1046  aXPLUS_bY(dat1.data(), Yptr, opts.a_arg, opts.b_arg,
1047  h1_slice_comps,h1->num_components(), out);
1048  break;
1049  case op_arg_aXPLUS_b:
1050  aXPLUS_b(dat1.data(),opts.a_arg,opts.b_arg,h1_slice_comps,out);
1051  break;
1052  case op_arg_XdivY:
1053  if( Ytype==SPACE_TIME || Ytype==SPACE )
1054  XdivY(dat1.data(), Yptr, h1_slice_comps, out);
1055  else
1056  XdivY(dat1.data(), Yptr, h1_slice_comps, h1->num_components(), out);
1057  break;
1058  case op_arg_XmulY:
1059  if( Ytype==SPACE_TIME || Ytype==SPACE )
1060  XmulY( dat1.data(), Yptr, h1_slice_comps, out);
1061  else
1062  XmulY( dat1.data(), Yptr, h1_slice_comps, h1->num_components(), out);
1063  break;
1064  case op_arg_minX_t:
1065  minXt( dat1.data(), h1->slice_sz(), h1->num_components(), tm, out );
1066  break;
1067  case op_arg_maxX_t:
1068  maxXt( dat1.data(), h1->slice_sz(), h1->num_components(), tm, out );
1069  break;
1070  case op_arg_avgX_t:
1071  avgXt( dat1.data(), h1->slice_sz(), h1->num_components(), tm, out );
1072  break;
1073  case op_arg_minX_i:
1074  minXi( dat1.data(), t_reduce.data(), t_reduce.size() );
1075  break;
1076  case op_arg_maxX_i:
1077  maxXi( dat1.data(), t_reduce.data(), t_reduce.size() );
1078  break;
1079  case op_arg_t_maxX:
1080  t_maxX( dat1.data(), big_scratch.data(), t_reduce.data(), t_reduce.size(), t );
1081  break;
1082  case op_arg_avgX_i:
1083  case op_arg_stdD_i:
1084  avgXi( dat1.data(), t_reduce, 1.0f);
1085  break;
1086  case op_arg_ravgX_ta :
1087  ravgX_ta( dat1.data(), h1_slice_comps, t, opts.a_arg, big_scratch, out );
1088  break;
1089  case op_arg_Xbox :
1090  Xbox( dat1.data(), h1->x(), h1->y(), out, opts.a_arg, h1->num_components() );
1091  break;
1092  case op_arg_sclip_ab :
1093  sclip_ab( dat1.data(), h1_slice_comps, opts.a_arg, opts.b_arg, big_scratch, t, h1->t(), out );
1094  break;
1095  case op_arg_diff_f :
1096  forward_diff( dat1.data(), h1_slice_comps, big_scratch, t, h1->inc_t(), out );
1097  break;
1098  case op_arg_diff_c :
1099  center_diff( dat1.data(), h1_slice_comps, big_scratch, t, h1->inc_t(), out );
1100  break;
1101  case op_arg_joinXY_i:
1102  if( Yigb )
1103  join_space( dat1.data(), h1_slice_comps, Yptr,
1104  h2->slice_sz()*h2->num_components(), out);
1105  else
1106  join_space( dat1.data(), h1_slice_comps, Y.data(), Y.size(), out);
1107  break;
1108  case op_arg_magX:
1109  magnitude( dat1.data(), h1_slice_comps, h1->num_components(), out );
1110  break;
1111  case op_arg_dotXY:
1112  dot( dat1.data(), Yptr, h1_slice_comps, h1->num_components(), out,
1113  Ytype==PT_TIME|| Ytype==PT );
1114  break;
1115  case op_arg_intX_i :
1116  case op_arg_intX_t :
1117  {
1118  int lastT = h1->t()-1;
1119  int lastSimp = lastT - (lastT%3);
1120  if( t<=lastSimp && lastSimp!=0 )
1121  simp38( dat1.data(), h1_slice_comps, t==0||t==lastSimp?-1:t, t_reduce.data(),
1122  h1->inc_t() );
1123  if( t>=lastSimp && lastSimp!=lastT )
1124  midpt( dat1.data(), h1_slice_comps, t==lastSimp||t==lastT, t_reduce.data(),
1125  h1->inc_t() );
1126  }
1127  if( opts.op_arg == op_arg_intX_t )
1128  fwrite( t_reduce.data(), sizeof(t_reduce[0]), t_reduce.size(), out );
1129  break;
1130  case op_arg_norm_i:
1131  minXi( dat1.data(), big_scratch.data(), h1_slice_comps );
1132  maxXi( dat1.data(), big_scratch.data()+h1_slice_comps, h1_slice_comps );
1133  break;
1134  case op_arg_maxP2P:
1135  maxP2P( dat1.data(), big_scratch.data(), t_reduce.data(), h1_slice_comps, t );
1136  break;
1137  case op_arg_NOP:
1138  break;
1139  }
1140  }
1141 
1142  // post process
1143  switch( opts.op_arg ) {
1144  case op_arg_avgX_i:
1145  for( int i=0; i<h1->slice_sz()*h1->num_components(); i++ )
1146  t_reduce[i] /= (data_t)(h1->t());
1147  break;
1148  case op_arg_stdD_i:
1149  {
1150  for( int i=0; i<h1->slice_sz()*h1->num_components(); i++ )
1151  t_reduce[i] /= (data_t)(h1->t());
1152  std::vector<data_t> mean;
1153  mean.assign(t_reduce.begin(), t_reduce.end());
1154  std::fill(t_reduce.begin(), t_reduce.end(), 0);
1155  for( int t=0; t<h1->t(); t++ ) {
1156  float tm = h1->org_t() + t*h1->inc_t();
1157  varXi( dat3[t].data(), mean.data(), t_reduce, 1.0f);
1158  }
1159  for( int i=0; i<h1->slice_sz()*h1->num_components(); i++ )
1160  t_reduce[i] = std::sqrt(t_reduce[i] / ((data_t)(h1->t() - 1)));
1161  break;
1162  }
1163  case op_arg_norm_i:
1164  normalize( h1, out, big_scratch.data(), big_scratch.data()+h1_slice_comps );
1165  break;
1166  case op_arg_t_maxX:
1167  for( int i=0; i<t_reduce.size(); i++ )
1168  t_reduce[i] = h1->org_t() + t_reduce[i]*h1->inc_t();
1169  break;
1170  }
1171 
1172  if( t_reduce.size() )
1173  output_scratch( t_reduce, out, true );
1174 
1175  return 0;
1176 }
constexpr T min(T a, T b)
Definition: ion_type.h:33
int type(void)
Definition: IGBheader.h:280
void read_vector(vector< data_t > &y, gzFile in, int nc)
Definition: IGBops.cc:742
int read_data(T *dp, size_t numt=1, char *buf=NULL)
Definition: IGBheader.h:431
void XdivY(data_t *x, data_t *y, int slsz, FILE *out)
Definition: IGBops.cc:334
void midpt(T *x, int n, bool end, data_t *reduce, float dt)
Definition: IGBops.cc:594
void normalize(IGBheader *h, FILE *out, T *min, T *max)
normalize pixels over all time
Definition: IGBops.cc:637
void aXPLUS_bY(data_t *x, data_t *y, data_t a, data_t b, int slsz, int nc, FILE *out)
Definition: IGBops.cc:274
void magnitude(data_t *x, int n, int nc, FILE *out)
determine magnitude at each point
Definition: IGBops.cc:529
#define IGB_FLOAT
Definition: IGBheader.h:60
void eval_expr(data_t *x, int n, char *expr, FILE *out, data_t *result)
Definition: IGBops.cc:172
constexpr T max(T a, T b)
Definition: ion_type.h:31
size_t slice_sz()
Definition: IGBheader.h:263
int go2slice(size_t s)
go to a slice
Definition: IGBheader.cc:1459
void center_diff(data_t *x, int slsz, vector< data_t > &odat, int t, float dt, FILE *out)
Definition: IGBops.cc:239
size_t t(void)
Definition: IGBheader.h:277
#define IGB_VEC3_f
Definition: IGBheader.h:71
vector< int * > ind
Definition: IGBops.cc:81
void aXPLUS_b(data_t *x, data_t a, data_t b, int slsz, FILE *out)
Definition: IGBops.cc:325
void interpolate(data_t *dat, vector< int *> &idx, vector< float *> &w, int n, FILE *out)
interpolate within points
Definition: IGBops.cc:616
Interpolate class.
Definition: IGBops.cc:76
#define data_t
Definition: IGBops.cc:121
Definition: IGBops.cc:70
void forward_diff(data_t *x, int slsz, vector< data_t > &old, int t, float dt, FILE *out)
Definition: IGBops.cc:261
void avgXi(data_t *x, vector< data_t > &sum, const data_t fac)
Definition: IGBops.cc:308
int Tbit
Definition: IGBops.cc:71
void maxXi(data_t *x, data_t *mx, size_t n)
Definition: IGBops.cc:423
void output_scratch(data_t *s, FILE *out, int n, bool bin)
Definition: IGBops.cc:694
void XPLUS_a(data_t *x, data_t a, int slsz, FILE *out)
Definition: IGBops.cc:458
const int MAX_NC
Definition: IGBops.cc:66
unsigned short num_components()
Definition: IGBheader.h:378
void joinXY(IGBheader *h1, IGBheader *h2, IGBheader *hout)
Definition: IGBops.cc:135
float org_t(void)
Definition: IGBheader.h:317
void XmulY(data_t *x, data_t *y, int slsz, FILE *out)
Definition: IGBops.cc:379
Definition: IGBops.cc:70
int read(bool quiet=false)
Definition: IGBheader.cc:757
int main(int argc, char *argv[])
Definition: IGBops.cc:771
int n
Definition: IGBops.cc:82
void ravgX_ta(data_t *dat1, int sz, int t, int a, vector< data_t > &tmp, FILE *out)
Definition: IGBops.cc:718
void maxXt(data_t *x, int slsz, int nc, float t, FILE *out)
Definition: IGBops.cc:445
std::string cmdline_str(int n, char *w[])
Definition: IGBops.cc:759
T sum(const vector< T > &vec)
Compute sum of a vector&#39;s entries.
Definition: SF_vector.h:340
void join_space(data_t *x, int nx, data_t *y, int ny, FILE *out)
Definition: IGBops.cc:164
Interp(char *ind, char *weight)
constructor
Definition: IGBops.cc:91
void varXi(data_t *x, data_t *mean, vector< data_t > &sum, const data_t fac)
Definition: IGBops.cc:316
bool overlap(IGBheader *first, IGBheader *second)
Definition: IGBops.cc:128
void inc_t(float a)
Definition: IGBheader.h:329
void maxP2P(data_t *dat, data_t *big, data_t *maxpp, int sz, int t)
find the maximum peak-to-peak value at each point
Definition: IGBops.cc:668
void simp38(T *x, int n, int t, data_t *runInt, float dt)
Simpson&#39;s 3/8 rule.
Definition: IGBops.cc:572
size_t x(void)
Definition: IGBheader.h:268
int IGB_DATA_T[5]
Definition: IGBops.cc:122
int size()
Definition: IGBops.cc:79
void t_maxX(data_t *x, data_t *mx, data_t *tm, size_t n, float t)
Definition: IGBops.cc:411
void dot(data_t *x, data_t *y, int n, int nc, FILE *out, bool pt)
perform dot product
Definition: IGBops.cc:550
void Xbox(data_t *dat, int x, int y, FILE *out, int d, int nc)
Definition: IGBops.cc:469
void avgXt(data_t *x, int slsz, int nc, float t, FILE *out)
Definition: IGBops.cc:293
void sclip_ab(data_t *dat, int sz, data_t upper, data_t lower, vector< data_t > &scr, int t, int tfinal, FILE *out)
Definition: IGBops.cc:493
void minXi(data_t *x, data_t *mn, size_t n)
Definition: IGBops.cc:434
size_t y(void)
Definition: IGBheader.h:271
Definition: IGBops.cc:70
void fileptr(gzFile f)
Definition: IGBheader.cc:329
#define IGB_VEC4_f
Definition: IGBheader.h:73
Data_file_t
Definition: IGBops.cc:70
#define IGB_COMPLEX
Definition: IGBheader.h:62
vector< float * > w
Definition: IGBops.cc:80
double mag(const Point &vect)
vector magnitude
Definition: SF_container.h:111
void minXt(data_t *x, int slsz, int nc, float t, FILE *out)
Definition: IGBops.cc:397
Definition: IGBops.cc:70