二维和三维平滑滤波（Gauss平滑、SU Smooth、 Dave Hale‘s method）

Gauss平滑计算量大，在大模型平滑时非常耗时，在此基础上，提供另外两种高效的平滑程序，分别来自于SU软件和Dave Hale的递归滤波平滑的修改版本。

/* 3 kinds smooth  */ 
#include <iostream>
#include <fstream>
#include <math.h>
#include "tool.h"
using namespace std;


#if 0
/***** Smooth from Hale*****/
reference to hale's home:
https://inside.mines.edu/~dhale/notebook.html
#endif
int smooth1Ei(bool zs,float m,int n,float *x,float *y);
int smooth2Ei(bool zs,float m,int n1,int n2,float **x,float **y);
int smooth3Ei(bool zs,float m,int n1,int n2,int n3,float ***x,float ***y);
int smooth3Ei1(bool zs,float m,int n1,int n2,int n3,float ***x,float ***y);




#if 0
/***** Smooth from SU*****/
SMOOTH3D - 3D grid velocity SMOOTHing by the damped least squares	

smooth3d <infile >outfile [parameters]				

Required Parameters:							
n1=		number of samples along 1st dimension			
n2=		number of samples along 2nd dimension			
n3=		number of samples along 3rd dimension			

Optional Parameters:							

Smoothing parameters (0 = no smoothing)				
	r1=0.0	operator length in 1st dimension			
	r2=0.0	operator length in 2nd dimension			
	r3=0.0	operator length in 3rd dimension			

	Sample intervals:							
	d1=1.0	1st dimension						
	d2=1.0	2nd dimension						
	d3=1.0	3rd dimension						
	iter=2	number of iteration used				
	time=0	which dimension the time axis is (0 = no time axis)	
	depth=1	which dimension the depth axis is (ignored when time>0)	
mu=1		the relative weight at maximum depth (or time)		
	verbose=0	=1 for printing minimum wavelengths			
	slowness=0	=1 smoothing on slowness; =0 smoothing on velocity	
	vminc=0	velocity values below it are cliped before smoothing	
	vmaxc=99999	velocity values above it are cliped before smoothing	

	Notes:								
	1. The larger the smoothing parameters, the smoother the output velocity.
	These parameters are lengths of smoothing operators in each dimension.
	2. iter controls the orders of derivatives to be smoothed in the output
	velocity. e.g., iter=2 means derivatives until 2nd order smoothed.	
	3. mu is the multipler of smoothing parameters at the bottom compared to
	those at the surface.						
	4. Minimum wavelengths of each dimension and the total may be printed	
	for the resulting output velocity is. To compute these parameters for
	the input velocity, use r1=r2=r3=0.					
	5. Smoothing directly on slowness works better to preserve traveltime.
	So the program optionally converts the input velocity into slowness	", 
	and smooths the slowness, then converts back into velocity.		

	Author:  CWP: Zhenyue Liu  March 1995

	Reference:
	Liu, Z., 1994,"A velocity smoothing technique based on damped least squares
	in Project Reveiew, May 10, 1994, Consortium Project on
	Seismic Inverse Methods for Complex Stuctures.
#endif
	void vsm3d(float ***v,int n3,int n2,int n1,int iter,int depth,float r3,float r2,float r1,float mu,int sl,float vmin,float vmax);
	void tripd2(float **d, float **e, float **b, int n, int m);


/*Gauss smooth*/
void smoothgauss3d(float ***fv,float ***fvf,int n1,int n2,int n3,int dx,int dy,int dz,int wx,int wy,int wz,float **elev);





int main(void)
{
	int nx=1573;
	int ny=498;
	int nz=75;
	float ***data=new float**[nx];
	float **sf=new float*[nx];
	float ***Edata=new float**[nx];
	//ifstream fin("V_YJZ.dat0");//../data/Record/shot_gather.dat_CorreT_loop0_
	//ifstream fin("MA01d_HQ3D_LP3_TSEQ_FB1_ALL_Rb1_Rb1_Srn_O7_3000_ini_extend.datSmall");//../data/Record/shot_gather.dat_CorreT_loop0_
	ifstream fin("dv.dat0");//../data/Record/shot_gather.dat_CorreT_loop0_
	for(int ix=0;ix<nx;ix++)
	{
		data[ix]=new float*[ny];
		sf[ix]=new float[ny];
		Edata[ix]=new float*[ny];
		for(int iy=0;iy<ny;iy++)
		{
			data[ix][iy]=new float [nz];
			Edata[ix][iy]=new float [nz];
			for(int iz=0;iz<nz;iz++)
			{ 
				data[ix][iy][iz]=0;	
				Edata[ix][iy][iz]=0;	
				fin.read((char*)&data[ix][iy][iz],sizeof(float));
				Edata[ix][iy][iz]=data[ix][iy][iz];
			}
		}

	}
	fin.close();
	for(int ix=0;ix<nx;ix++)
	{
		for(int iy=0;iy<ny;iy++)
		{
			for(int iz=0;iz<nz;iz++)
			{
				if(data[ix][iy][iz]>340) 
				{
					sf[ix][iy]=iz;
					break;
				}
			}
		}
	}
	/*	for(int ix=0;ix<nx;ix++)
		for(int iy=0;iy<41;iy++)
		for(int iz=0;iz<nz;iz++)
		Edata[ix][iy][iz] = Edata[ix][41][iz];
		for(int ix=0;ix<nx;ix++)
		for(int iy=301;iy<ny;iy++)
		for(int iz=0;iz<nz;iz++)
		Edata[ix][iy][iz] = Edata[ix][300][iz];
		*/
	//smooth2Ei(true,30,ny,nx,data1,data1);
	//smooth3Ei(true,50,nz,ny,nx,Edata,Edata);
	//smooth3Ei(false,30,nz,ny,nx,Edata,Edata);
	//smoothgauss3d(data,Edata,nx,ny,nz,20,20,10,4,4,2,sf);

	vsm3d(Edata,nx,ny,nz,2,1,0.5*30*30,0.5*30*30,0.5*10*10,1,0,-9999,3400);
	//void vsm3d(float ***v,int n3,int n2,int n1,int iter,int depth,float r3,float r2,float r1,float mu,int sl,float vmin,float vmax);
	ofstream fout("Smooth.dat");
	for(int ix=0;ix<nx;ix++)
	{
		for(int iy=0;iy<ny;iy++)
		{
			for(int iz=0;iz<nz;iz++)
			{
				fout.write((char*)&Edata[ix][iy][iz],sizeof(float));
			}
		}
	}
	fout.close();
	for(int ix=0;ix<nx;ix++)
	{
		for(int iy=0;iy<ny;iy++)
		{

			delete [] data[ix][iy];
			delete [] Edata[ix][iy];
		}
		delete [] data[ix];
		delete [] Edata[ix];
		delete [] sf[ix];
	}
	delete [] data;
	delete [] Edata;
	delete [] sf;

	return 0;
}



int smooth1Ei(bool zs,float m,int n,float *x,float *y)
{
	float sigma = sqrt(m*(m+1)/3.0);
	float ss = sigma*sigma;
	float a = (1.0+ss-sqrt(1.0+2.0*ss))/ss;
	float b=1.0-a;
	float sx = zs?1.0:b;
	float sy = a;
	float yi = 0.0;
	yi = y[0] = sx*x[0];
	for(int i=1;i<n-1;++i)
		y[i] = yi = a*yi+b*x[i];
	sx /= 1.0+a;
	sy /=1.0+a;
	y[n-1] = yi = sy*yi+sx*x[n-1];
	for(int i=n-2;i>=0;--i)
		y[i] = yi = a*yi+b*y[i];
}//end of smooth1Ei()
//Smooth along 2nd of a 2D array for input boundary conditions(Hale)
int smooth2Ei(bool zs,float m,int n1,int n2,float **x,float **y)
{
	float sigma = sqrt(m*(m+1)/3.0);
	float ss = sigma*sigma;
	float a = (1.0+ss-sqrt(1.0+2.0*ss))/ss;
	float b = 1-a;
	float sx = zs?1.0:b;
	float sy = a;

//	for(int i=0;i<n2;i++)
//		smooth1Ei(zs,n1/10,n1,x[i],x[i]);

	for(int i1=0;i1<n1;++i1) 
		y[0][i1] = sx*x[0][i1];

	for(int i2=1;i2<n2-1;++i2)
		for(int i1=0;i1<n1;++i1)
			y[i2][i1] = a*y[i2-1][i1]+b*x[i2][i1];

	sx /= 1+a;
	sy /= 1+a;

	for(int i1=0;i1<n1;++i1)
		y[n2-1][i1] = sy*y[n2-2][i1]+sx*x[n2-1][i1];

	for(int i2=n2-2;i2>=0;--i2)
		for(int i1=0;i1<n1;++i1)
			y[i2][i1] = a*y[i2+1][i1]+b*y[i2][i1];
}//end of smooth2Ei()
int smooth3Ei1(bool zs,float m,int n1,int n2,int n3,float ***x,float ***y)
{
	for(int i3=0;i3<n3;++i3)
	{
		smooth2Ei(true,30,n1,n2,x[i3],x[i3]);
		smooth2Ei(true,30,n1,n2,y[i3],y[i3]);
	}
}//end of smooth3Ei()
//Smooth along 3nd of a 3D array for input boundary conditions(Extended for Hale 2D)
int smooth3Ei(bool zs,float m,int n1,int n2,int n3,float ***x,float ***y)
{
	float sigma = sqrt(m*(m+1)/3.0);
	float ss = sigma*sigma;
	float a = (1.0+ss-sqrt(1.0+2.0*ss))/ss;
	float b = 1-a;
	float sx = zs?1.0:b;
	float sy = a;

	for(int i=0;i<n3;i++){
		for(int j=0;j<n2;j++){
			smooth1Ei(zs,10,n1,x[i][j],x[i][j]);
		}
	}
	for(int i=0;i<n3;i++)
		smooth2Ei(zs,30,n1,n2,x[i],x[i]);


	for(int i2=0;i2<n2;++i2) 
		for(int i1=0;i1<n1;++i1) 
			y[0][i2][i1] = sx*x[0][i2][i1];

	for(int i3=1;i3<n3-1;++i3)
		for(int i2=0;i2<n2;++i2)
			for(int i1=0;i1<n1;++i1)
				y[i3][i2][i1] = a*y[i3-1][i2][i1]+b*x[i3][i2][i1];

	sx /= 1+a;
	sy /= 1+a;

	for(int i2=0;i2<n2;++i2)
		for(int i1=0;i1<n1;++i1)
			y[n3-1][i2][i1] = sy*y[n3-2][i2][i1]+sx*x[n3-1][i2][i1];

	for(int i3=n3-2;i3>=0;--i3)
		for(int i2=0;i2<n2;++i2)
			for(int i1=0;i1<n1;++i1)
				y[i3][i2][i1] = a*y[i3+1][i2][i1]+b*y[i3][i2][i1];

	//add by zjm
	//	 for(int i1=0;i1<n1;++i1) 
	//		for(int i3=0;i3<n3;++i3) 
	//		{
	//			y[i3][0][i1] = y[i3][1][i1];
	//			y[i3][n2-1][i1] = y[i3][n2-2][i1];
	//		}

}//end of smooth3Ei()

void vsm3d(float ***v,int n3,int n2,int n1,int iter,int depth,float r3,float r2,float r1,float mu,int sl,float vmin,float vmax)
	/***************************************************************************
	 * Smooth 3d-velocity.  
	 * *************************************************************************/

{
	int  i2, i1, i3, i;		
	float **d=NULL, **e=NULL, **f=NULL, *w, ww=1.0;

	/*	compute the weight function */
	//w = alloc1float(n1+n2+n3-2);
	w = alloc_float_1d(n1+n2+n3-2);
	if(depth==1){
		mu = (mu*mu-1.0)/(n1*n1);
		for(i1=0; i1<n1; ++i1) w[i1] = 1.0/(1+i1*i1*mu);
	}
	if(depth==2){
		mu = (mu*mu-1.0)/(n2*n2);
		for(i2=0; i2<n2; ++i2) w[i2] = 1.0/(1+i2*i2*mu);
	}
	if(depth==3){
		mu = (mu*mu-1.0)/(n3*n3);
		for(i3=0; i3<n3; ++i3) w[i3] = 1.0/(1+i3*i3*mu);
	}

	/*	scale  smoothing parameters according to the iteration
	 *	number	*/
	if(iter==1) {
		r1 /= 3.39*3.39;
		r2 /= 3.39*3.39;
		r3 /= 3.39*3.39;
	} else if(iter==2){
		r1 /= 5.19*5.19;
		r2 /= 5.19*5.19;
		r3 /= 5.19*5.19;
	} else {
		r1 /= 6.60*6.60;
		r2 /= 6.60*6.60;
		r3 /= 6.60*6.60;
	}


	/*  clip velocity  */
	for(i3=0; i3<n3; ++i3) 
		for(i2=0; i2<n2; ++i2)
			for(i1=0; i1<n1; ++i1){
				if(v[i3][i2][i1] >vmax) v[i3][i2][i1] = vmax;
				if(v[i3][i2][i1] <vmin) v[i3][i2][i1] = vmin;
			}

	if(sl) {
		/*  smoothing on slowness  */
		for(i3=0; i3<n3; ++i3) 
			for(i2=0; i2<n2; ++i2)
				for(i1=0; i1<n1; ++i1)
					v[i3][i2][i1] = 1.0/v[i3][i2][i1];
	}


	if(r2>0.) {

		/*	smoothing velocity in the second
		 *	direction */

		/* allocate space */
		//d = alloc2float(n1,n2);
		//e = alloc2float(n1,n2);
		//f = alloc2float(n1,n2);
		d = alloc_float_2d(n2,n1);
		e = alloc_float_2d(n2,n1);
		f = alloc_float_2d(n2,n1);


		for(i3=0; i3<n3; ++i3){
			if(depth==3) ww = w[i3];
			for(i2=0; i2<n2-1; ++i2){
				if(depth==2) ww = w[i2+1];
				for(i1=0; i1<n1; ++i1){
					if(depth==1) ww = w[i1];
					d[i2][i1] = ww+r2*2.0;
					e[i2][i1] = -r2;
					f[i2][i1] = ww*v[i3][i2+1][i1];
				}
			}
			for(i1=0; i1<n1; ++i1){
				d[n2-2][i1] -= r2;
				f[0][i1] += r2*v[i3][0][i1];
			}
			tripd2(d,e,f,n2-1,n1);

			for(i=1; i<iter; ++i) {
				for(i2=0; i2<n2-1; ++i2){
					if(depth==2) ww = w[i2+1];
					for(i1=0; i1<n1; ++i1){
						if(depth==1) ww = w[i1];
						d[i2][i1] = ww+r2*2.0;
						e[i2][i1] = -r2;
						f[i2][i1] *= ww;
					}
				}
				for(i1=0; i1<n1; ++i1){
					d[n2-2][i1] -= r2;
					f[0][i1] += r2*v[i3][0][i1];
				}
				tripd2(d,e,f,n2-1,n1);
			}

			for(i2=0; i2<n2-1; ++i2)
				for(i1=0; i1<n1; ++i1)
					v[i3][i2+1][i1] = f[i2][i1];
		}
	}

	if(r3>0.) {
		/*	smooth velocity in  the third
		 *	direction */

		/* allocate space */
		//d = alloc2float(n1,n3);
		//e = alloc2float(n1,n3);
		//f = alloc2float(n1,n3); 
		d = alloc_float_2d(n3,n1);
		e = alloc_float_2d(n3,n1);
		f = alloc_float_2d(n3,n1);

		for(i2=0; i2<n2; ++i2){
			if(depth==2) ww = w[i2];
			for(i3=0; i3<n3-1; ++i3){
				if(depth==3) ww = w[i3+1];
				for(i1=0; i1<n1; ++i1){
					if(depth==1) ww = w[i1];
					d[i3][i1] = ww+2.*r3;
					e[i3][i1] = -r3;
					f[i3][i1] = ww*v[i3+1][i2][i1];
				}
			}
			for(i1=0; i1<n1; ++i1){
				d[n3-2][i1] -= r3;
				f[0][i1] += r3*v[0][i2][i1];
			}
			tripd2(d,e,f,n3-1,n1);

			for(i=1; i<iter; ++i){
				for(i3=0; i3<n3-1; ++i3){
					if(depth==3) ww = w[i3+1];
					for(i1=0; i1<n1; ++i1){
						if(depth==1) ww = w[i1];
						d[i3][i1] = ww+2.*r3;
						e[i3][i1] = -r3;
						f[i3][i1] *= ww;
					}
				}
				for(i1=0; i1<n1; ++i1){
					d[n3-2][i1] -= r3;
					f[0][i1] += r3*v[0][i2][i1];
				}
				tripd2(d,e,f,n3-1,n1);
			}

			for(i3=0; i3<n3-1; ++i3)
				for(i1=0; i1<n1; ++i1)
					v[i3+1][i2][i1] = f[i3][i1];
		}
	}

	if(r1>0.) {
		/*	smooth velocity in  the
		 *	first direction */

		/* allocate space */
		//d = alloc2float(1,n1);
		//e = alloc2float(1,n1);
		//f = alloc2float(1,n1);
		d = alloc_float_2d(n1,1);
		e = alloc_float_2d(n1,1);
		f = alloc_float_2d(n1,1);

		for(i3=0; i3<n3; ++i3){
			if(depth==3) ww = w[i3];
			for(i2=0; i2<n2; ++i2){
				if(depth==2) ww = w[i2];
				for(i1=0; i1<n1-1; ++i1){
					if(depth==1) ww = w[i1+1];
					d[i1][0] = ww+r1*2.0;
					e[i1][0] = -r1;
					f[i1][0] = ww*v[i3][i2][i1+1];
				}
				d[n1-2][0] -= r1;
				f[0][0] += r1*v[i3][i2][0];
				tripd2(d,e,f,n1-1,1);

				for(i=1; i<iter; ++i) {
					for(i1=0; i1<n1-1; ++i1){
						if(depth==1) ww = w[i1+1];
						d[i1][0] = ww+r1*2.0;
						e[i1][0] = -r1;
						f[i1][0] *= ww;
					}
					d[n1-2][0] -= r1;
					f[0][0] += r1*v[i3][i2][0];
					tripd2(d,e,f,n1-1,1);
				}

				for(i1=0; i1<n1-1; ++i1)
					v[i3][i2][i1+1] = f[i1][0];
			}
		}
	}

	if(sl) {
		for(i3=0; i3<n3; ++i3) 
			for(i2=0; i2<n2; ++i2)
				for(i1=0; i1<n1; ++i1)
					v[i3][i2][i1] = 1.0/v[i3][i2][i1];
	}

	//free1float(w);
	free_1d(w);
	if(r1>0. || r2>0. || r3>0.) {
		//free2float(d);
		//free2float(e);
		//free2float(f);
		free_2d(d,n1);
		free_2d(e,n1);
		free_2d(f,n1);
	}
}
void tripd2(float **d, float **e, float **b, int n, int m)
	/*****************************************************************************
	 * Given m n-by-n symmetric, tridiagonal, positive definite matri2 A's and m
	 * n-vector b's, the following algorithm overwrites b with the solution to Ax = b.
	 * The first dimension of arrays is independent of the algorithm. 
	 *
	 *   d() the diagonal of A 
	 *     e() the superdiagonal of A
	 *     *****************************************************************************/
{
	int k, i; 
	float temp;

	/* decomposition */
	for(k=1; k<n; ++k) 
		for(i=0; i<m; ++i){ 
			temp = e[k-1][i];
			e[k-1][i] = temp/d[k-1][i];
			d[k][i] -= temp*e[k-1][i];
		}

	/* substitution	*/
	for(k=1; k<n; ++k) 
		for(i=0; i<m; ++i) 
			b[k][i] -= e[k-1][i]*b[k-1][i];

	for(i=0; i<m; ++i) 
		b[n-1][i] /=d[n-1][i];

	for(k=n-1; k>0; --k) 
		for(i=0; i<m; ++i) 
			b[k-1][i] = b[k-1][i]/d[k-1][i]-e[k-1][i]*b[k][i]; 

}



void smoothgauss3d(float ***fv,float ***fvf,int n1,int n2,int n3,int dx,int dy,int dz,int wx,int wy,int wz,float **elev)
{
	int i,j,m,i2,j2,m2;
	int lx,ly,lz,k1,k2,k3;
	float betatot;

	lx=3*wx;
	ly=3*wy;
	lz=3*wz;

	float *betax,*betay,*betaz;
	betax =new float[2*lx+1];
	betay =new float[2*ly+1];
	betaz =new float[2*lz+1];


	for(i=0;i<2*lx+1;i++)
		betax[i]=exp(-((i-lx)*dx)*(i-lx)*dx/(wx*dx*wx*dx));
	for(i=0;i<2*ly+1;i++)
		betay[i]=exp(-((i-ly)*dy)*(i-ly)*dy/(wy*dy*wy*dy));
	for(i=0;i<2*lz+1;i++)
		betaz[i]=exp(-((i-lz)*dz)*(i-lz)*dz/(wz*dz*wz*dz));

	for(i=0;i<n1;i++)
		for(j=0;j<n2;j++)
			for(m=0;m<n3;m++)
			{
				if(m<elev[i][j])
					fvf[i][j][m]=fv[i][j][m];
				else{
					fvf[i][j][m]=0;
					betatot=0;
					for(i2=0;i2<2*lx+1;i2++)
					{
						k1=i+(i2-lx);
						if(k1>=0 && k1<n1)
						{
							for(j2=0;j2<2*ly+1;j2++)
							{
								k2=j+(j2-ly);
								if(k2>=0 && k2<n2)
								{
									for(m2=0;m2<2*lz+1;m2++)
									{	
										k3=m+(m2-lz);
										if(k3>=elev[i][j] && k3<n3)
										{
											fvf[i][j][m]+=betax[i2]*betay[j2]*betaz[m2]*fv[k1][k2][k3];
											betatot+=betax[i2]*betay[j2]*betaz[m2];
										}
									}
								}
							}
						}
					}
					fvf[i][j][m]=fvf[i][j][m]/betatot;
				}
			}

}