supickamp.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* SUPICKAMP: $Revision: 1.9 $ ; $Date: 2006/10/31 22:20:54 $	*/

#include "su.h"
#include "segy.h"


/*********************** self documentation **********************/
char *sdoc[] = {
" 	   								",
" SUPICKAMP - pick amplitudes within user defined and resampled window	",
" 	   								",
"   supickamp <stdin >stdout d2=  [optional parameters]			",
" 									",
" Required parameters:							",
" d2=		   sampling interval for slow dimension			",
"			(required if key-parameter not specified)	",
" Optional parameters:							",
" key= 			Key header word specifying trace offset 	",
" 			(alternatively, specify d2,x2beg)		",
" 									",
" x_above=		array of lateral position values   		",
" 			(upper window corner)				",
" t_above=		array of time values   				",
" 			(upper window corner)				",
"									",
"  ... or input via files:						",
" t_xabove=		file containing time and lateral position values",
" 			(upper window corner)				",
" t_xbelow=		file containing time and lateral position values",
" 			(lower window corner)				",
" wl=		   	window width if t_xbelow is not specified	",
"			(No windowing if not specified)		 	",
" 									",
" dt_resamp=dt	  resampling interval within pick window	  	",
"			(dt has to come from trace headers)		",
" tmin=0		minimum time in input trace			",
" x2beg=0		first lateral position				",
" format=ascii 		write ascii data to stdout			",
"			 =binary for binary floats to stdout		",
" verbose=1 		writes complete  pick information into outpar   ",
"			=2 writes complete pick information into outpar	",
"			   in tab-delimited column format		",
" outpar=/dev/tty	output parameter file; contains output		",
"					from verbose			",
" arg1=max		output (first dimension) to stdout		",
" arg2=i2		output (second dimension) to stdout		",
"			(see notes for other options)			",
" Notes: 								",
"									",
" Window can be defined using						",
" (1)   vectors x_above, t_above, [wl]					",
" (2)   file  t_xabove, [wl]	or					",
" (3)   files t_xabove, t_xbelow					",
"									",
" files t_xabove, t_xbelow can be generated using xwigb's picking	",
" algorithm. The lateral positions have to be monotonically increasing  ",
" or decreasing for both vector and file input.				",
" verbose=1 or 2 writes min, max, abs[max], energy and associated times	",
" tmin,tmax,tabs to outpar, together with global values. verbose=0	",
" only outputs global values.						",
" Acceptable arg-parameters for lateral positions are			",
" (1) x2   (2) i2 = trace number					",
"									",
" If key=keyword is set, then the values of x2 are taken from the header",
" field represented by the keyword (for example key=offset)		",
" Type	sukeyword -o   to see the complete list of SU keywords.		",
"									",
NULL};

/*
 * Credits:
 *
 *	CWP: Andreas Rueger July 06, 1996
 *	MTU: David Forel,   Jan. 26, 2005  Add verbose=2 option
 */
/**************** end self doc ***********************************/

/* Locally defined macros */
#define diprint(expr) printf(#expr " = %i\n",expr)
#define dfprint(expr) printf(#expr " = %f\n",expr)
#define ddprint(expr) printf(#expr " = %g\n",expr)

segy tr;

/* Pick parameterss */
typedef struct PickStruct {
	float max,min,abs;
	float tmax,tmin,tabs;
	float energy;
	int imax,imin,iabs,ienergy;
} Pick;

/* Prototype of functions used internally */
int tracepick(Pick *pick1, float *t, float *val, int nt);
int initpick(Pick *pick1);
int globalpick(Pick *pick1, Pick *gpick);
int get_mode(char *arg, int *mode);
int getformat(char *format, int *modef);
int getarg(float *arg, int mode, Pick *pick, float x2,int trnum);


int
main(int argc, char **argv)
{
	char *arg1,*arg2;	/* control parameters		*/
	char *key;		/* pointer SEGY keyword		*/
	char *type=NULL;	/* ... its type			*/
	int index=0;		/* ... its index		*/
	Value val;		/* ... its value		*/
	float fval;		/* ... its value cast to float	*/
	
   	float *x_above=NULL;	/* pointer to x values above	*/
	float *x_below=NULL;	/* pointer to x values below	*/
	float *t_above=NULL;	/* pointer to t values above	*/
	float *t_below=NULL;	/* pointer to t values above	*/

	float wl;		/* window width			*/
	float x2beg;		/* first  x2 value		*/		
	float d2;		/* x2 sampling interval		*/
	float t1;		/* time values			*/
	float x1;		/* x values			*/
	float dt_resamp,tmin,dt,tmin_in;
	
	float *t_resamp=NULL;	/* resampled time		*/
	float *val_resamp=NULL;	/* resampled amplitude		*/
	float tmax;		/* maximum time			*/
	
	int nx,i,trnum;	
	int nt,nts,modef;	
	int verbose,mode1,mode2 ;

	char *outpar=NULL;	/* name of file holding output parfile  */
	FILE *outparfp=NULL;	/* ... its file pointer			*/

	char *format=NULL;	/* format (ascii/binary) of output to stdout */

	/* Names of input files */
	cwp_String afile="";	
	FILE *afilep=NULL;		
	cwp_String bfile="";	
	FILE *bfilep=NULL;		

	/* Booleans to keep track of open files */
	cwp_Bool afile_set=cwp_false;
	cwp_Bool bfile_set=cwp_false;

	Pick *pick1,*gpick;
	
	/* Initialize */
	initargs(argc, argv);
	requestdoc(1);

	/* Allocate space for picks */
	pick1=(Pick*)emalloc(1*sizeof(Pick));
	gpick=(Pick*)emalloc(1*sizeof(Pick));


	/* Get parameters */
	if (!getparfloat("d2", &d2)) d2=0; 
	if (getparstring("key", &key))	{
	   d2=0;
	} else if(d2==0) {
	   err(" Specify d2");
	}

	if (!getparint("verbose", &verbose))	verbose = 1;
	if (!getparstring("arg1", &arg1))	arg1 = "max";
	if (!getparstring("arg2", &arg2))	arg2 = "x2";

	if (!getparfloat("wl", &wl))	wl = 0.;
	if (!getparfloat("dt_resamp", &dt_resamp))   dt_resamp = 0.;
	if (!getparfloat("x2beg", &x2beg))	x2beg = 0.;
	if (!getparfloat("tmin", &tmin_in))	tmin_in = 0.;

	if (!getparstring("format", &format)) format="ascii";

	if (!getparstring("outpar", &outpar))  outpar = "/dev/tty" ;
	outparfp = efopen(outpar, "w");

	/* column output */
	if(verbose == 2)
	fprintf(outparfp,
		"trace\tx\tt_above\tt_below\t"
		"trace #\tx2\t"
		"max\t tmax\t"
		"min\t tmin\t"
		"abs\t tabs\t"
		"energy in pick window\n") ;

	/* get key type and index */
	if(d2==0.) {
	   type = hdtype(key);
	   index = getindex(key);
	}


	if (getparstring("t_xabove",&afile)) {

		/* check if file exists */
		if ((afilep=fopen(afile,"r")) == NULL)
		fprintf(stderr,"can't open t_xabove file \n");

		/* Quick runthrough */
		nx=0;
   
		while (fscanf(afilep,"%f %f \n",&t1,&x1 ) != EOF) ++nx;

		rewind(afilep);

		x_above=ealloc1float(nx);
		t_above=ealloc1float(nx);  

		/* Read above_file */
		nx=0;
			while (fscanf(afilep,"%f %f \n",&t1,&x1 ) != EOF) {
				x_above[nx]=x1;
				t_above[nx]=t1;

				/* Test output 
				fprintf(outparfp,"%f %f \n",t1,x1); */
		
				++nx;
			}
		

		} else {

			nx=nt=0;
		
			nx = countparval("x_above");
			nt = countparval("t_above");

		if (nx != nt)
			err("lengths of x_above, t_above must be the same");

		x_above = ealloc1float(nx);  getparfloat("x_above", x_above);
		t_above = ealloc1float(nx);  getparfloat("t_above", t_above);

		/* Test output 

		for(i=0;i<nx;i++)
			fprintf(outparfp,"%f %f \n",t_above[i],x_above[i]); */

	}

	/* Lower limit picks of time-window */
	if (getparstring("t_xbelow",&bfile) ) {

		if (afilep == NULL)
		err(" Specify t_xabove -file");
		
		/* check if file exists */
		if ((bfilep=fopen(bfile,"r")) == NULL)
		fprintf(stderr,"can't open t_xbelow-file \n");


		/* Quick runthrough */
		nx=0;
   
		while (fscanf(bfilep,"%f %f \n",&t1,&x1 ) != EOF)
		{
		   nx++;
		}

		rewind(bfilep);

		x_below=ealloc1float(nx);
		t_below=ealloc1float(nx);  

		/* Read below_file */
		nx=0;
			while (fscanf(bfilep,"%f %f \n",&t1,&x1 ) != EOF)
		{
		x_below[nx]=x1;
		t_below[nx]=t1;

		/* Test output 
		fprintf(outparfp,"%f %f \n",t1,x1); */
		
		nx++;
		}

	 }
	
	/* Specify window width */
	if(bfilep!= NULL && afilep != NULL) {
		afile_set = cwp_true;
		bfile_set = cwp_true;
	   warn("window determined by above_ and below_ files ");
	} else if(afilep != NULL && wl != 0) {
		afile_set = cwp_true;
	   warn("window determined by above_file and wl");
	} else if(afilep != NULL && wl == 0 ){
		afile_set = cwp_true;
	   warn("window determined by above_file and trace end");
	} else if(afilep == NULL && wl != 0 && nx!=0) { 
		afile_set = cwp_false;
	   warn("window determined by above_vector and wl");
	} else if(afilep == NULL && wl == 0 && nx!=0) {
		afile_set = cwp_false;
	   warn("window determined by above_vector and trace length");
	} else if(afilep == NULL && wl == 0) {
		afile_set = cwp_false;
	   warn("window determined by full trace length");
	} else {
	   err("window not properly specified");
	}


	/* get mode for stdout information */
	if(! get_mode(arg1,&mode1) )
	   err(" Unknown arg1 ");

	if(! get_mode(arg2,&mode2) )
	   err(" Unknown arg2");

	/* get format for stdout information */
	if(! getformat(format, &modef) )
	   err(" Unknown format"); 

	/* Get info from first trace */
	if (!gettr(&tr)) err("can't read first trace");
	if (!tr.dt) err("dt header field must be set");
	if (tr.delrt)
	   tmin_in = tr.delrt/1000.0;

	nt	= (int) tr.ns;
	dt   = ((double) tr.dt)/1000000.0;

	if(dt_resamp == 0)
	   dt_resamp=dt;
	
	/* allocate space for max search window */
	nts = (int) (nt-1)*dt/dt_resamp +1.;

	/* aux output 
	diprint(nts);
	diprint(nt);