openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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