C++ sf_floatwrite函数代码示例

void pt2dwrite1(sf_file F, 
		pt2d   *v, 
		size_t  n1, 
		int     k)
/*< output point2d 1-D vector >*/
{
    int i1;
    float **w;
    w=sf_floatalloc2(k,n1);

    for( i1=0; i1<(int)n1; i1++) {
	;        w[i1][0] = v[i1].x; 
	;        w[i1][1] = v[i1].z; 
	if(k==3) w[i1][2] = v[i1].v;
    }
    sf_floatwrite(w[0],k*n1,F);

    free( *w); free(w);
}

int main (int argc, char* argv[])
{
    bool scale;
    int n1,n2, i1,i2, n;
    float d1, *dat, *der;
    sf_file in=NULL, out=NULL;

    sf_init (argc,argv);
    in = sf_input("in");
    out = sf_output("out");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    n2 = sf_leftsize(in,1);

    dat = sf_floatalloc(n1);
    der = sf_floatalloc(n1);

    if (!sf_getint("order",&n)) n=6;
    /* filter order */

    if (!sf_getbool("scale",&scale) || 
	!sf_histfloat(in,"d1",&d1)) scale=false;
    /*(scale=n if scale by 1/dx )*/

    sf_deriv_init(n1, n, 0.);

    for (i2=0; i2 < n2; i2++) {
	sf_floatread(dat,n1,in);
	sf_deriv(dat,der);

	if (scale) {
	    for (i1=0; i1 < n1; i1++) {
		der[i1] /= d1;
	    }
	}

	sf_floatwrite(der,n1,out);
    }


    exit(0);
}

int main(int argc, char* argv[])
{
    bool adj, inv;
    int nn,n1,n2,i1,i2;
    float rho, *pp=NULL;
    sf_file in=NULL, out=NULL;

    sf_init(argc,argv);
    in = sf_input("in");
    out = sf_output("out");

    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    n2 = sf_leftsize(in,1);

    if (!sf_getbool("adj",&adj)) adj=false;
    /* If y, apply adjoint */
    if (!sf_getbool("inv",&inv)) inv=false;
    /* If y, do differentiation instead of integration */
    if (!sf_getfloat("rho",&rho)) rho = 1.-1./n1;
    /* Leaky integration constant */

    if (inv) {
	sf_warning("%s half-order differentiation",adj? "anticausal":"causal");
    } else {
	sf_warning("%s half-order integration",adj? "anticausal":"causal");
    }

    nn = 2*kiss_fft_next_fast_size((n1+1)/2);
    pp = sf_floatalloc(nn);

    sf_halfint_init (inv, nn, rho);

    for (i2=0; i2 < n2; i2++) {
	sf_floatread (pp,n1,in);
	for (i1=n1; i1 < nn; i1++) {
	    pp[i1]=0.;
	}
	sf_halfint (adj, pp);
	sf_floatwrite (pp,n1,out);
    }

    exit(0);
}

int main(int argc, char* argv[])
{
    int nt, ns, nx, ix, is, it;
    float *trace=NULL, *picks=NULL, **semb=NULL, s0, ds, s;
    sf_file in=NULL, out=NULL, pick=NULL;

    sf_init(argc,argv);
    in = sf_input("in");
    pick = sf_input("pick");
    out = sf_output("out");

    if (!sf_histint(in,"n1",&nt)) sf_error("No n1= in input");

    if (!sf_histint(in,"n2",&ns)) sf_error("No n2= in input");
    if (!sf_histfloat(in,"o2",&s0)) sf_error("No o2= in input");
    if (!sf_histfloat(in,"d2",&ds)) sf_error("No d2= in input");

    nx = sf_unshiftdim(in,out,2);

    semb = sf_floatalloc2(nt,ns);
    picks = sf_floatalloc (nt);
    trace = sf_floatalloc (nt);

    for (ix=0; ix < nx; ix++) {
	sf_floatread (semb[0],nt*ns,in);
	sf_floatread (picks,nt,pick);

	for (it=0; it < nt; it++) {
	    s = (picks[it] - s0)/ds;
	    is = floorf(s); s -= is;
	    if (is >= 0 && is < ns-1) {
		trace[it] = s * semb[is+1][it] + (1.-s) * semb[is][it];
	    } else {
		trace[it] = 0.;
	    }
	}

	sf_floatwrite (trace,nt,out);
    }

    exit (0);
}

int main(int argc, char* argv[])
{
    int  i1, i2, n1,n2,n3,n12;
    float  s=0, *pp1, *pp2;
    sf_file inp1, inp2,  dif;

    sf_init(argc,argv);
    inp1 = sf_input("in");
    inp2 = sf_input("match");
    dif = sf_output("out");
    
    sf_putfloat(dif,"o2",0);
    sf_putfloat(dif,"o1",0);
    sf_putint(dif,"n1",1);  
    sf_putint(dif,"n2",1);    
	
    if (!sf_histint(inp1,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(inp1,"n2",&n2)) sf_error("No n2= in input");
    if (!sf_histint(inp1,"n3",&n3)) n3=1;

    n2 = sf_leftsize(inp1,1); /* left dimensions after the first one */
    n12 = n1*n2;
    
    if (n3!=1) {sf_putint(dif,"n3",1); sf_putint(dif,"d3",1); }
    
    pp1 = sf_floatalloc(n12);
    pp2 = sf_floatalloc(n12);
    sf_floatread(pp1,n12,inp1);
    sf_floatread(pp2,n12,inp2);

    for (i2=0; i2 < n2; i2++) {
	for (i1=0; i1 < n1; i1++) {
	    s = s + pow((pp1[i2*n1+i1]-pp2[i2*n1+i1]),2);
	}
    }
    
    sf_warning("The difference is %f", s );
    sf_floatwrite(&s,1,dif);
    exit(0);
}

int main(int argc, char* argv[])
{
    int nx, nz, ix, iz;
    float* trace;
    float dx, dz, x[2];
    sf_file mod;

    sf_init (argc,argv);
    mod = sf_output("out");

    if (!sf_getint("nx",&nx)) nx=400; 
    /* horizontal dimension */
    if (!sf_getint("nz",&nz)) nz=800; 
    /* vertical dimension */

    dx = 1./(nx-1);
    dz = 2./(nz-1);

    sf_putint   (mod,"n1",nz); 
    sf_putfloat (mod,"d1",dz); 
    sf_putfloat (mod,"o1",0.);
    sf_putint   (mod,"n2",nx); 
    sf_putfloat (mod,"d2",dx); 
    sf_putfloat (mod,"o2",0.);
    sf_setformat (mod,"native_float");

    trace = sf_floatalloc(nz);

    for (ix = 0; ix < nx; ix++) {
	x[1] = ix*dx;
	for (iz = 0; iz < nz; iz++) {
	    x[0] = iz*dz;
	    trace[iz] = 1./sqrtf(symes_vel(NULL,x));
	}
	sf_floatwrite(trace,nz,mod);
    }

    exit (0);
}

void pt2dwrite2(sf_file F, 
		pt2d  **v, 
		size_t  n1, 
		size_t  n2,
		int     k)
/*< output point2d 2-D vector >*/
{
    int i1,i2;
    float ***w;
    w=sf_floatalloc3(k,n1,n2);

    for( i2=0; i2<(int)n2; i2++) {
	for( i1=0; i1<(int)n1; i1++) {
	    ;        w[i2][i1][0] = v[i2][i1].x; 
	    ;        w[i2][i1][1] = v[i2][i1].z; 
	    if(k==3) w[i2][i1][2] = v[i2][i1].v;
	}
    }
    sf_floatwrite(w[0][0],k*n1*n2,F);
    
    free(**w); free(*w); free(w);
}

int main(int argc, char* argv[])
{
    int n1, n2, i2;
    float *pp, *qq;
    bool square, adj;
    sf_file in, out;

    sf_init(argc,argv);

    in = sf_input("in");
    out = sf_output("out");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    n2 = sf_leftsize(in,1);

    if (!sf_getbool("square",&square)) square=false;
    /* if y, use gradient squared */

    if (!sf_getbool("adj",&adj)) adj=false;
    /* adjoint flag */

    pp = sf_floatalloc(n1);
    qq = sf_floatalloc(n1);

    for (i2=0; i2 < n2; i2++) {
	sf_floatread(pp,n1,in);
	if (square) {
	    sf_grad2 (n1,pp,qq);
	} else if (adj) {
	    sf_igrad1_lop (true,false,n1,n1,qq,pp);
	} else {
	    sf_igrad1_lop (false,false,n1,n1,pp,qq);
	}
	sf_floatwrite(qq,n1,out);
    }

    exit(0);
}

int main(int argc,char *argv[])
{
    int i,j,n1,n2,n3;
    float o1,d1,o2,d2,max;
    float *orig,*output;
    sf_file in,out;
    
    sf_init (argc, argv); 
    in = sf_input("in");
    out = sf_output("out");
    
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histfloat(in,"d1",&d1)) sf_error("No d1= in input");
    if (!sf_histfloat(in,"o1",&o1)) sf_error("No o1= in input");
    if (!sf_histint(in,"n2",&n2)) n2=1;
    if (!sf_histfloat(in,"d2",&d2)) d2=1;
    if (!sf_histfloat(in,"o2",&o2)) o2=1;
    n3=sf_leftsize(in,2);
    
    orig = sf_floatalloc(n1*n2);
    output = sf_floatalloc(n1*n2);
    for(j = 0; j < n3; j++){
	
	sf_floatread(orig,n1*n2,in);
	max=0;
	for(i=0;i<n1*n2;i++){
	    if(orig[i]>max) max=orig[i];
	}
	
	for(i=0;i<n1*n2;i++){
	    output[i]=max-orig[i];
	}
	
	sf_floatwrite(output,n1*n2,out);
    }
    exit(0);
    
}

int main(int argc, char* argv[])
{ 
    int dim, m[SF_MAX_DIM], nd, i, niter, rect[SF_MAX_DIM];
    float **inp, *rat;
    char key[6];
    sf_file in, out;

    sf_init (argc, argv);
    in = sf_input("in");
    out = sf_output("out");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    dim = sf_filedims(in,m);

    nd=1;
    for (i=0; i < dim; i++) {
	nd *= m[i];
	snprintf(key,6,"rect%d",i+1);
	if (!sf_getint(key,rect+i)) rect[i]=1;
    }

    if (!sf_getint("niter",&niter)) niter=100;
    /* number of iterations */
    
    inp = sf_floatalloc2 (m[0],nd/m[0]);
    rat = sf_floatalloc (nd);

    envcorr_init(dim,m,rect,niter);
    
    sf_floatread(inp[0],nd,in);

    envcorr(inp,rat);

    sf_floatwrite(rat,nd,out);

    exit (0);
}

int main(int argc, char** argv)
{
    int n1, n2, i1, i2;
    float* vint=NULL;
    float o1, d1, v0, alpha, t;
    sf_file in=NULL, out=NULL;

    sf_init (argc, argv);
    in = sf_input("in");
    out = sf_output("out");

    if (!sf_histint (in,"n1",&n1)) sf_error ("No n1= in input\n");
    if (!sf_histint (in,"n2",&n2)) sf_error ("No n2= in input\n");

    if (!sf_histfloat (in,"o1",&o1)) o1=0.;
    if (!sf_histfloat (in,"d1",&d1)) sf_error ("No d1= in input\n");

    if (!sf_getfloat ("v0",&v0)) v0 = 1.5;
    /* initial velocity */
    if (!sf_getfloat ("alpha",&alpha)) alpha = 0.5;
    /* velocity gradient */

    vint = sf_floatalloc(n1);
    for (i1=0; i1 < n1; i1++) {
	t = alpha*(o1+i1*d1);
	if (t > 0.) {
	    vint[i1] = v0 * sqrtf((expf(t)-1.)/t);
	} else {
	    vint[i1] = v0;
	}
    }
    for (i2=0; i2 < n2; i2++) {
	sf_floatwrite(vint,n1,out);
    }

    exit (0);
}

int main (int argc, char *argv[])
{
    int n1, n2, n3, i2, i3, order;
    float *trace;
    sf_file in, out;

    sf_init(argc,argv);
    in = sf_input("in");
    out = sf_output("out");

    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (!sf_histint(in,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2)) sf_error("No n2= in input");
    n3 = sf_leftsize(in,2);
    
    if (!sf_getint("order",&order)) order=4;
    /* Accuracy order */

    trace = sf_floatalloc (n1);

    pick0_init (n1, n2, order);

    for (i3=0; i3 < n3; i3++) {
	for (i2=0; i2 < n2; i2++) {
	    sf_floatread(trace,n1,in);
	    pick0_set (i2, trace);
	}

	for (i2=0; i2 < n2; i2++) {
	    pick0_delta (i2, trace);
	    sf_floatwrite(trace,n1,out);
	}
    }
    

    exit (0);
}

int main(int argc, char *argv[])
{
    int i, m, n, k, n1a, n1b, n2a, n2b, n3a, n3b; /*n1 is trace length, n2 is the number of traces, n3 is the number of 3th axis*/
    sf_file in, mat, out;
    float *a, *b, *c;
    char *type;
    sf_init(argc,argv);

    in =  sf_input("in");    
    /*  in is the input matrix A, with m rows and n columns */
    mat= sf_input("mat");
    /* mat is the input matrix B, with n rows and k columns */
    out = sf_output("out");
    /* out is the output matrix C, with m rows and k columns */

    if (NULL == (type=sf_getstring("type"))) type="mul";
    /* [mul, add, sub, dotmul] operation type, the default is mul  */


    if (!sf_histint(in,"n1",&n1a)) sf_error("No n1= in input");
    if (!sf_histint(in,"n2",&n2a)) sf_error("No n2= in input");
    n3a = sf_leftsize(in,2);

    if (!sf_histint(mat,"n1",&n1b)) sf_error("No n1= in input");
    if (!sf_histint(mat,"n2",&n2b)) sf_error("No n2= in input");
    n3b = sf_leftsize(mat,2);
    
    if (type[0]=='m')
	{
	if (n1a!=n2b) sf_error("n1 of A must be equal to n2 of B for multiplication");
	m=n2a;n=n1a;k=n1b;
	}

    if (type[0]!='m')
	{
	if (n1a!=n1b || n2a!=n2b) sf_error(" n1 of A must be equal to n1 of B and n2 of A must be equal to n2 of B");
	m=n2a;n=n1a;k=n;	
	}
    if (n3a != n3b) sf_error("n3 must be equal");
    
    sf_putint(out, "n1", k);   
    sf_putint(out, "n2", m);
    sf_putint(out, "n3", n3a);

    a = sf_floatalloc(m*n);
    b = sf_floatalloc(n*k);
    c = sf_floatalloc(m*k);

    for(i=0;i<n3a;i++)  {
	sf_floatread(a, m*n, in);

	if(type[0]!='m')
	sf_floatread(b, m*n, mat);
	else
	sf_floatread(b,n*k,mat);

	switch(type[0])
	{
		case 'm': mul(a, b, m, n, k, c); break;
		case 'a': add(a, b, m, n, c); break;
		case 's': sub(a, b, m, n, c); break;
		case 'd': dotmul(a, b, m, n, c);
	}
	sf_floatwrite(c, m*k, out);
    }

    exit(0);
}

int main (int argc, char *argv[])
{
    bool verb;
    int n1,n2,n3, n12, n23, ref2, ref3, i2,i3,i1, ud, lr, order;
    float eps, ***dat, ***p, ***q, **p2, **q2, *trace;
    sf_file dip, out, seed, cost;

    sf_init(argc,argv);
    dip = sf_input("in");
    out = sf_output("out");
    seed = sf_input("seed");
    cost = sf_input("cost");

    if (!sf_histint(dip,"n1",&n1)) sf_error("No n1= in input");
    if (!sf_histint(dip,"n2",&n2)) sf_error("No n2= in input");
    if (!sf_histint(dip,"n3",&n3)) sf_error("No n3= in input");
    n23 = n2*n3;
    n12 = n1*n23;

    if (!sf_getbool("verb",&verb)) verb=false;
    if (!sf_getfloat("eps",&eps)) eps=0.01;
    /* regularization */

    if (!sf_getint("ref2",&ref2)) ref2=0;
    if (!sf_getint("ref3",&ref3)) ref3=0;
    /* reference trace */

    sf_putint(out,"n4",1);

    p = sf_floatalloc3(n1,n2,n3);
    q = sf_floatalloc3(n1,n2,n3);
    dat = sf_floatalloc3(n1,n2,n3);

    sf_floatread(p[0][0],n12,dip);
    sf_floatread(q[0][0],n12,dip);

    p2 = sf_floatalloc2(n2,n3);
    q2 = sf_floatalloc2(n2,n3);

    sf_floatread(p2[0],n23,cost);
    sf_floatread(q2[0],n23,cost);

    dijkstra_init(n2,n3,p2,q2);
    dijkstra_source(ref2,ref3);
    sf_floatread(dat[ref3][ref2],n1,seed);
    if (verb) sf_warning("%d %d",ref2,ref3);

    if (!sf_getint("order",&order)) order=1;
    /* accuracy order */

    predict_init(n1,n2, eps*eps, order, 1, false);
    
    while (dijskstra_step(&i2,&i3,&ud,&lr)) {
	if (verb) sf_warning("%d %d",i2,i3);

	trace = dat[i3][i2];
	for (i1=0; i1 < n1; i1++) {
	    trace[i1] = dat[i3-lr][i2-ud][i1];
	} 
	if (ud > 0) {
	    predict_step(false,true,trace,p[i3][i2-ud]);
	} else if (ud < 0) {
	    predict_step(false,false,trace,p[i3][i2]);
	}

	if (lr > 0) {
	    predict_step(false,true,trace,q[i3-lr][i2]);
	} else if (lr < 0) {
	    predict_step(false,false,trace,q[i3][i2]);
	}
    }
	
    sf_floatwrite(dat[0][0],n12,out);

    exit (0);
}

//.........这里部分代码省略.........
    sf_putfloat(out,"o2",o2); 
    sf_putfloat(out,"o3",0.0); 

    nxb = nx + nxl + nxr;
    nzb = nz + nbt + nbb;
    
	sf_putint(out,"n1",nzb);
    sf_putint(out,"n2",nxb);
    sf_putint(out,"n3",nsnap);


    wav    =  sf_floatalloc(nt);
    sf_floatread(wav,nt,source);

    old    =  sf_floatalloc2(nzb,nxb);
    cur    =  sf_floatalloc2(nzb,nxb);
    aa     =  sf_floatalloc3(nzb,nxb,3);
    

    bd2_init(nx,nz,nxl,nxr,nbt,nbb,cxl,cxr,czt,czb);

    /*input & extend velocity model*/
    v = sf_floatalloc2(nzb,nxb);
    vtmp = sf_floatalloc2(nz,nx);
    sf_floatread(vtmp[0],nx*nz,vel);
	v = extmodel(vtmp, nz, nx, nbt);

    sigma = sf_floatalloc2(nzb,nxb);
    sf_floatread(vtmp[0],nx*nz,fsigma);
	sigma = extmodel(vtmp, nz, nx, nbt);

    delta = sf_floatalloc2(nzb,nxb);
    sf_floatread(vtmp[0],nx*nz,fdelta);
	delta = extmodel(vtmp, nz, nx, nbt);
    
    seta = sf_floatalloc2(nzb,nxb);
    sf_floatread(vtmp[0],nx*nz,fseta);
	seta = extmodel(vtmp, nz, nx, nbt);

    v0 =0.0;
    for (ix=0; ix < nxb; ix++) {
        for (iz=0; iz < nzb; iz++) {
            v0 += v[ix][iz]*v[ix][iz];
         }
    }

    v0 = sqrtf(v0/(nxb*nzb));

	fprintf(stderr, "v0=%f\n\n", v0);


    for (ix=0; ix < nxb; ix++) {
        for (iz=0; iz < nzb; iz++) {
            cur[ix][iz] = 0.0;
            old[ix][iz] = 0.0; 
        }
    }

    /* propagation in time */
    pamstep2_init(nzb,nxb,dz,dx,opt);

    for (it=0; it < nt; it++) {
		fprintf(stderr, "\b\b\b\b\b%d", it);

        pamstep2(old, cur, nzb, nxb, dz, dx, v0, v, sigma, delta, seta, dt); 
        old[isx+nxl][isz+nbt] += wav[it];
     
		bd2_decay(old); 
        bd2_decay(cur); 
        tmp = old;
        old = cur;
        cur = tmp;

		if (it%nsnap==0)
		sf_floatwrite(cur[0], nxb*nzb, out);

    }

    pamstep2_close();
    bd2_close();

    free(**aa);
    free(*aa);
    free(aa);
    free(*v);     
    free(*sigma);     
    free(*delta);     
    free(*seta);     
	free(*vtmp);
    free(*cur);     
    free(*old);     
    free(v);     
    free(sigma);     
    free(delta);     
    free(seta);     
	free(vtmp);
    free(cur);     
    free(old);     
    exit(0); 
}