smoothint2.c

su 的源代码库

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

/* SMO0THINT2: $Revision: 1.8 $ ; $Date: 2003/08/19 21:24:44 $	*/

#include "par.h"

/*********************** self documentation ******************************/
char *sdoc[] = {
" 									",
" SMOOTHINT2 --- SMOOTH non-uniformly sampled INTerfaces, via the damped",
"  		least-squares technique					",
" 									",
"  smoothint2 <input ninf= >output [optional parameters]		",
" 									",
" Required Parameters:							",
" <input                 file containing original interfaces		",
" >output                file containing smoothed interfaces	 	",
" 									",
" Optional Parameters:							",
" ninf=5                number of interfaces  				",
" r=100			smoothing parameter 				",
" npmax=101		maximum number of points in interfaces		",
" 									",
" Notes:								",
" The input file is an ASCII file. Each interface is represented by pairs",
" (non-uniform sampling) of x and z values, with one pair of values on	",
" each line, separated by spaces or tabs. Each interface is separated with",
" an entry with a large negative z value for example: 1.0     -9999.	",
" There is no entry for the surface. The surface is assumed to be flat  ",
" with z=0.								",
" This is similar to a CSHOT model file without a surface entry and	",
" without comments.							",
" 									",
" The smoothing method is analogous to a moving window averaging process",
" (but not the same) with the parameter \"r\" being analogous to the \"width\"",
" of the window. Thus, the size of \"r\" must be chosen to by compatible",
" with the scale (wavelengths) of the variations of the interfaces in the",
" model being smoothed.							",
" 									",
" Example using the test data set generated by unif2: 			",
" unif2 tfile=tfilename							",
" Compare the unsmoothed interface model:				",
" unif2 < tfilename method=interpolation_method |			",
"	 			psimage n1=100 n2=100 d1=10 d2=10 | ...	",
" To the smoothed interface model:					",
" smoothint2 r=100 < tfilename | unif2 method=interpolation_method |	", 
"	psimage n1=100 n2=100 d1=10 d2=10 | ...				",
"   									",
NULL};

/*
 * Credits:
 *  CWP: Zhenyue Liu, Jan 1994
 * Reference:
 *    Liu, Zhenyue, 1994, Velocity smoothing: theory and implementation, 
 *    Project Review, 1994, Consortium Project on Seismic Inverse Methods
 *    for Complex Stuctures (in review)
 */

/**************** end self doc *******************************************/

void smooth_1(float *x, float *z, float r, int n);

int
main(int argc, char **argv)
{
	int i,j,ninf,npmax,npt;
	float r;
	float *xint,*zint;
	
	/* hook up getpar */
	initargs(argc, argv);
	requestdoc(0);
	
	/* get parameters */
	if (!getparint("ninf",&ninf))	ninf=5;
	if (!getparint("npmax",&npmax)) npmax=101;
	if (!getparfloat("r",&r)) r=0;
	
	xint =  alloc1float(npmax);
	zint =  alloc1float(npmax);
	
	/* Input x z pairs specifying interfaces and velocity values. */
	for(i=0; i<ninf; ++i) {
		j= -1;
		do {
			j +=1;
			if(j>=npmax) 
		err("The point number on the %ith interface exceeds %i!",
				i+1,npmax);
			scanf("%f%f\n", &xint[j], &zint[j]);
		} while(zint[j]>-9999);
			npt = j;

		smooth_1(xint,zint,r,npt);
		for(j=0; j<npt; ++j)
 			printf("\t%g\t%g\n",xint[j],zint[j]);
 		printf("\t%g\t%g\n",1.,-99999.);
	}

	return(CWP_Exit());
}	       

void smooth_1(float *x, float *z, float r, int n)
/**************************************************************************
smooth_1  --- damped least squares smoothing of nonuniformly sampled 1d array
**************************************************************************
Input:
	x     array of horizontal coordinates
	z     array of vertical coordinate
	r     smoothing operator length
	n     number of points input 
Output:
	z     array of smoothed vertical coordinates
**************************************************************************
Reference:
    Liu, Zhenyue, 1994, Velocity smoothing: theory and implementation, 
    Project Review, 1994, Consortium Project on Seismic Inverse Methods
    for Complex Stuctures (in review).
**************************************************************************
Author: CWP: Zhenyue Liu, 1994
**************************************************************************/
{
	int  ir, i;		
	float *d, *e, *o;

	o = alloc1float(n+1);
	d = alloc1float(n);
	e = alloc1float(n);

	for(i=0; i<n-1; ++i){
		if(x[i+1]==x[i])
 			err("x-spacing of data values must not be zero!\n");
		o[i] = 1.0/(x[i+1]-x[i]);
		o[i] = o[i]*o[i];
	}
	o[n-1] = 0.;
		
	/* scale smoothing parameter */
	r = 0.5*r*r;
	r /= 5.19*5.19;

	for(ir=0; ir<2; ++ir) {
		for(i=1; i<n; ++i){
			d[i] = 1.0+r*(o[i]+o[i-1]);
			e[i-1] = -r*o[i-1];
		}
		d[0] = 1.0+r*o[0];
         	tripd(d,e,z,n);
	}

	free1float(d);
	free1float(e);
	free1float(o);
}