sugain.c

su 的源代码库

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

/* SUGAIN: $Revision: 1.51 $ ; $Date: 2006/11/07 21:16:33 $               */

#include "su.h"
#include "segy.h"
#include "header.h"
#include <signal.h>

/*********************** self documentation *****************************/
char *sdoc[] = {
"                                                                       ",
" SUGAIN - apply various types of gain to display traces                ",
"                                                                       ",
" sugain <stdin >stdout [optional parameters]                           ",
"                                                                       ",
" Required parameters:                                                  ",
"       none (no-op)                                                    ",
"                                                                       ",
" Optional parameters:                                                  ",
"       panel=0         =1 gain whole data set (vs. trace by trace)	",
"       tpow=0.0        multiply data by t^tpow                         ",
"       epow=0.0        multiply data by exp(epow*t)                    ",
"       gpow=1.0        take signed gpowth power of scaled data         ",
"       agc=0           flag; 1 = do automatic gain control             ",
"       gagc=0          flag; 1 = ... with gaussian taper               ",
"       wagc=0.5        agc window in seconds (use if agc=1 or gagc=1)  ",
"       trap=none       zero any value whose magnitude exceeds trapval  ",
"       clip=none       clip any value whose magnitude exceeds clipval  ",
"       qclip=1.0       clip by quantile on absolute values on trace    ",
"       qbal=0          flag; 1 = balance traces by qclip and scale     ",
"       pbal=0          flag; 1 = bal traces by dividing by rms value   ",
"       mbal=0          flag; 1 = bal traces by subtracting the mean    ",
"       maxbal=0        flag; 1 = balance traces by subtracting the max ",
"       scale=1.0       multiply data by overall scale factor           ",
"       norm=1.0        divide data by overall scale factor             ",
"       bias=0.0        bias data by adding an overall bias value	",
"       jon=0           flag; 1 means tpow=2, gpow=.5, qclip=.95        ",
"       verbose=0       verbose = 1 echoes info				",
"       mark=0          apply gain only to traces with tr.mark=0	",
"                       =1 apply gain only to traces with tr.mark!=0    ",
"                                                                       ",
"                                                                       ",
" 	tmpdir=	 if non-empty, use the value as a directory path	",
"		 prefix for storing temporary files; else if the	",
"	         the CWP_TMPDIR environment variable is set use		",
"	         its value for the path; else use tmpfile()		",
"                                                                       ",
" Operation order:                                                      ",
"                                                                       ",
" out(t) = scale * BAL{CLIP[AGC{[t^tpow * exp(epow * t) * ( in(t)-bias )]^gpow}]}",
"                                                                       ",
" Notes:                                                                ",
"	The jon flag selects the parameter choices discussed in		",
"	Claerbout's Imaging the Earth, pp 233-236.			",
"									",
"	Extremely large/small values may be lost during agc. Windowing  ",
"	these off and applying a scale in a preliminary pass through	",
"	sugain may help.						",
"									",
"	Sugain only applies gain to traces with tr.mark=0. Use sushw,	",
"	suchw, suedit, or suxedit to mark traces you do not want gained.",
"	See the selfdocs of sushw, suchw, suedit, and suxedit for more	",
"	information about setting header fields. Use \"sukeyword mark\" ",
"	for more information about the mark header field.		",
"									",
NULL};

/* Credits:
 *      SEP: Jon Claerbout
 *      CWP: Jack K. Cohen, Brian Sumner, Dave Hale
 *
 * Note: Have assumed tr.deltr >= 0 in tpow routine.
 *
 * Technical Reference:
 *      Jon's second book, pages 233-236.
 *
 * Trace header fields accessed: ns, dt, delrt, mark
 */
/**************** end self doc *******************************************/

/* subroutine prototypes */
void gain(float *data, float tpow, float epow, float gpow,
	  int agc, int gagc, int qbal, int pbal, int mbal, float scale, float bias,
	  register float trap, register float clip, float qclip, int iwagc,
	  register float tmin, register float dt, int nt, int maxbal);
void do_tpow(float *data, float tpow, register float tmin, register float dt,
	     int nt);
void do_epow(float *data, float epow, register float tmin, register float dt,
	     int nt);
void do_trap(float *data, register float trap, register int nt);
void do_clip(float *data, register float clip, register int nt);
void do_qclip(float *data, float qclip, int nt);
void do_qbal(float *data, float qclip, int nt);
void do_agc(float *data, int iwagc, int nt);
void do_gagc(float *data, int iwagc, int nt);
float quant(float *a, int k, int n);
static void closefiles(void);

#define TPOW     0.0
#define EPOW     0.0
#define GPOW     1.0
#define TRAP     0.0
#define CLIP     0.0
#define QCLIP    1.0
#define SCALE    1.0
#define BIAS     0.0
#define WAGC     0.5

/* Globals (so can trap signal) defining temporary disk files */
char tracefile[BUFSIZ];	/* filename for the file of traces	*/
char headerfile[BUFSIZ];/* filename for the file of headers	*/
FILE *tracefp;		/* fp for trace storage file		*/
FILE *headerfp;		/* fp for header storage file		*/

segy tr;

int
main(int argc, char **argv)
{
        int verbose;	/* flag for echoing info			*/
	int jon;        /* flag to get Claerbout values                 */
	int agc;        /* agc flag                                     */
        int gagc;       /* gaussian agc flag                            */
        int pbal;       /* power balance flag                           */
        int qbal;       /* quantile balance flag                        */
        int mbal=0;     /* mean balance flag                            */
        float tpow;     /* exponent of t                                */
        float epow;     /* deattenutation coefficient                   */
        float gpow;     /* dynamic compression power                    */
        float trap;     /* zero any larger value magnitude than trapval */
        float clip;     /* clip any larger value magnitude than clipval */
        float qclip;    /* clip at qth quantile (100qth percentile)     */
        float scale;    /* overall scale factor                         */
        float norm;     /* reciprocal of scale factor                   */
        float bias=0.0; /* overall bias  value                          */
        float wagc;     /* size of agc window in seconds                */
        int iwagc=0;    /* ... half window in samples                   */
        int nt;         /* number of samples on trace                   */
        float tmin;     /* delay recording time in secs                 */
        float dt;       /* sample rate in secs                          */
	float *data;	/* the data					*/
	int panel;	/* gain trace by trace or whole data set?	*/

        int maxbal=0;   /* max balance flag                             */

        int mark;	/* mark flag					*/


	char *tmpdir;		/* directory path for tmp files		*/
	cwp_Bool istmpdir=cwp_false;/* true for user given path		*/


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

	/* Look for user-supplied tmpdir */
	if (!getparstring("tmpdir",&tmpdir) &&
	    !(tmpdir = getenv("CWP_TMPDIR"))) tmpdir="";
	if (!STREQ(tmpdir, "") && access(tmpdir, WRITE_OK))
		err("you can't write in %s (or it doesn't exist)", tmpdir);


	/* Get nt from first trace */
        if (!gettr(&tr)) err("can't get first trace");
        nt   = (int) tr.ns;
	dt = ((double) tr.dt)/1000000.0;   /* microsecs to secs */
	tmin = tr.delrt/1000.0;		   /* millisecs to secs */
        if (!dt) getparfloat("dt", &dt);
        if (!dt) MUSTGETPARFLOAT("dt", &dt);


        /* Get parameters */
        if (!getparfloat ("tpow" , &tpow))      tpow     = TPOW;
        if (!getparfloat ("epow" , &epow))      epow     = EPOW;
        if (!getparfloat ("gpow" , &gpow))      gpow     = GPOW;
        if (!getparfloat ("trap" , &trap))      trap     = TRAP;
        if (!getparfloat ("clip" , &clip))      clip     = CLIP;
        if (!getparfloat ("qclip", &qclip))     qclip    = QCLIP;
        if (!getparfloat ("scale", &scale))     scale    = SCALE;
        if (!getparfloat ("norm",  &norm)) 	norm     = 1.0;
        if (!getparfloat ("bias",  &bias))      bias     = BIAS;
        if (!getparfloat ("wagc" , &wagc))      wagc     = WAGC;
        if (!getparint   ("agc"  , &agc))       agc      = 0;
        if (!getparint   ("gagc" , &gagc))      gagc     = 0;
        if (!getparint   ("pbal" , &pbal))      pbal     = 0;
        if (!getparint   ("qbal" , &qbal))      qbal     = 0;
        if (!getparint   ("mbal" , &mbal))      mbal     = 0;
	if (!getparint   ("panel", &panel))     panel    = 0;
        if (!getparint   ("jon"  , &jon))       jon      = 0;
	if (!getparint("verbose", &verbose))	verbose  = 0;
        if (!getparint   ("maxbal" , &maxbal))      maxbal     = 0;
        if (!getparint   ("mark" , &mark))      mark     = 0;
	

        /* Data validation */
        if (trap < 0.0) err("trap = %f, must be positive", trap);
        if (clip < 0.0) err("clip = %f, must be positive", clip);
        if (qclip < 0.0 || qclip > 1.0) 
                err("qclip = %f, must be between 0 and 1", qclip);
        if (agc || gagc) {
                iwagc = NINT(wagc/dt);
                if (iwagc < 1) err("wagc=%g must be positive", wagc);
                if (iwagc > nt) err("wagc=%g too long for trace", wagc);
                iwagc >>= 1;  /* windows are symmetric, so work with half */
        }
        if (jon) { 
                tpow  = 2.0;
                gpow  = 0.5;
                qclip = 0.95;
        }

	if (norm==0.0) {
		err("parameter norm may not be zero!");
	} else {
		if (verbose && scale!=1.0)
			warn("value of scale overridden by norm!");
		scale = 1.0/norm;
	}

        /* Main loop over traces */
	if (!panel) { /* trace by trace */
		data = ealloc1float(nt);
		do {
			memcpy((void *)data, (const void *) tr.data, nt*FSIZE);

			if (!(tr.mark || mark) ) {
				gain(data, tpow, epow, gpow, agc, gagc, qbal,
				     pbal, mbal, scale, bias, trap, clip, qclip,
				     iwagc, tmin, dt, nt, maxbal);
			} else if ( (mark) && (tr.mark) ) {
				gain(data, tpow, epow, gpow, agc, gagc, qbal,
				     pbal, mbal, scale, bias, trap, clip, qclip,
				     iwagc, tmin, dt, nt, maxbal);
			}

			memcpy((void *)tr.data, (const void *) data, nt*FSIZE);
			puttr(&tr);

		} while(gettr(&tr));
	} else { /* do whole data set at once */
		int itr, ntr = 0;
		
		/* Store traces, headers in tempfiles while getting a count */
		if (STREQ(tmpdir,"")) {
			tracefp = etmpfile();
			headerfp = etmpfile();
			if (verbose) warn("using tmpfile() call");
		} else { /* user-supplied tmpdir */
		char directory[BUFSIZ];
		strcpy(directory, tmpdir);
		strcpy(tracefile, temporary_filename(directory));
		strcpy(headerfile, temporary_filename(directory));
		/* Handle user interrupts */
		signal(SIGINT, (void (*) (int)) closefiles);
		signal(SIGQUIT, (void (*) (int)) closefiles);
		signal(SIGHUP,  (void (*) (int)) closefiles);
		signal(SIGTERM, (void (*) (int)) closefiles);
		tracefp = efopen(tracefile, "w+");
		headerfp = efopen(headerfile, "w+");
      		istmpdir=cwp_true;		
		if (verbose) warn("putting temporary files in %s", directory);
	}
		do {
			++ntr;
			efwrite(&tr, HDRBYTES, 1, headerfp);
			efwrite(tr.data, FSIZE, nt, tracefp);
		} while (gettr(&tr));
		erewind(tracefp);
		erewind(headerfp);
		data = ealloc1float(nt*ntr);
		
		/* Load traces into data and close tmpfile */
		efread(data, FSIZE, nt*ntr, tracefp);
		efclose(tracefp);
		if (istmpdir) eremove(tracefile);
	
		gain(data, tpow, epow, gpow, agc, gagc, qbal,
		     pbal, mbal, scale, bias, trap, clip, qclip,
		     iwagc, tmin, dt, nt*ntr, maxbal);

		for (itr = 0; itr < ntr; itr++) {
			memcpy((void *) tr.data, (const void *) (data+itr*nt),
			       nt*FSIZE);
			efread(&tr, 1, HDRBYTES, headerfp);
			puttr(&tr);
		}
		efclose(headerfp);
		if (istmpdir) eremove(headerfile);
	}
	
	free1(data);
	
        return(CWP_Exit());
}


/* Multiply by t^tpow */
void do_tpow(
	float *data,		/* the data			*/
        float tpow,             /* multiply data by t^tpow      */
        register float tmin,    /* first time on record         */
        register float dt,      /* sampling rate in seconds     */
        int nt                  /* number of samples            */
)
{
        static cwp_Bool first = cwp_true;   /* first entry flag     */
        static float *tpowfac;          /* tpow values          */
        register int i;                 /* counter              */

        if (first) { /* first entry, set up array of tpow factors */
                tpowfac = ealloc1float(nt);

		/* protect against negative tpow */
		tpowfac[0] = (tmin == 0.0) ? 1.0 : pow(tmin, tpow);
                for (i = 1; i < nt; ++i)
                        tpowfac[i] = pow(tmin + i*dt, tpow);

                first = cwp_false;
        } /* end first entry */

        for (i = 0; i < nt; ++i)  data[i] *= tpowfac[i];
}


/* Exponential deattenuation  with deattenuation factor epow */
void do_epow(
	float *data,		/* the data			*/
        float epow,             /* coefficient of t in exponent */
        register float tmin,    /* first time on record         */
        register float dt,      /* sampling rate in seconds     */
        int nt                  /* number of samples            */
)
{
        register int i;                 /* counter              */
        static cwp_Bool first = cwp_true;   /* first entry flag     */
        static float *epowfac;          /* exponent stretchs    */

        if (first) {
                epowfac = ealloc1float(nt);
                for (i = 0; i < nt; i++) 
                        epowfac[i] = exp(epow * (tmin + i * dt));

                first = cwp_false;
        }

        for (i = 0; i < nt; ++i)  data[i] *= epowfac[i];
}


/* Zero out outliers */
void do_trap(
	float *data,		/* the data			*/
        register float trap,    /* zero if magnitude > trap     */
        register int nt         /* number of samples            */
)
{
        register float *dataptr = data;

        while (nt--) {
                if (ABS(*dataptr) > trap) *dataptr = 0.0;
                dataptr++;
        }
}


/* Hard clip outliers */
void do_clip(
	float *data,		/* the data				*/
        register float clip,    /* hard clip if magnitude > clip        */
        register int nt         /* number of samples                    */
)
{
        register float *dataptr = data;
	register float mclip = -clip;

        while (nt--) {
                if (*dataptr > clip) {