C++ sf_floatalloc3函数代码示例

void warp3_init(int n1_in, float o1, float d1,
		int n2_in, float o2, float d2,
		int n3_in, float o3, float d3 /* output grid */,
		int nt, int ny_in, int nx_in  /* input  grid */,
		float eps                     /* regularization */)
/*< initialize >*/
{
    n1 = n1_in;
    n2 = n2_in;
    n3 = n3_in;

    ny = ny_in;
    nx = nx_in;

    map1 = stretch4_init (n1, o1, d1, nt, eps);

    trace1 = sf_floatalloc(n1);
    trace2 = sf_floatalloc2(n2,n3);

    warp2_init(n2,o2,d2,
	       n3,o3,d3,
	       ny,nx,eps);
    
    str2   = sf_floatalloc3(ny,nx,n1);
    str3   = sf_floatalloc3(ny,nx,n1);
    slice1 = sf_floatalloc3(ny,nx,n1);
}

int main (int argc, char *argv[]) 
{    
    int m, k, l, seed;
    sf_file  bshuffle,ashuffle;
    sf_axis ax,at,ap,av;
    int nx,nt,np,nv,iteration,*a1, *a2;
    float ***bsh, ***bsh2;
    sf_init(argc,argv);

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

    if (!sf_getint("iteration",&iteration)) iteration=1;
    if (!sf_getint("seed",&seed)) seed=2012;
    at=sf_iaxa( bshuffle,1); nt=sf_n(at);
    ax=sf_iaxa( bshuffle,3); nx=sf_n(ax);
    ap=sf_iaxa( bshuffle,2); np=sf_n(ap);
    av=sf_iaxa( bshuffle,4); nv=sf_n(av);
    bsh=sf_floatalloc3(nt,np,nx);
    bsh2=sf_floatalloc3(nt,np,nx);
    a1=sf_intalloc(np);
    a2=sf_intalloc(np);
    sf_warning("ntpx=%d",nt*np*nx);

    srand(seed);


    for (m=0; m<np; m++) {
	a1[m]=rand();
    }

    bubble(a1, a2, np);


    for (k=0; k<nv; k++) {
	sf_floatread(bsh[0][0],nt*np*nx,bshuffle);
	for(l=0; l<nx; l++) {
	    for (m=0; m<np; m++) {
		memcpy(bsh2[l][m], bsh[l][a2[m]], nt*sizeof(float));
	    }
	}
	sf_floatwrite(bsh2[0][0],nt*np*nx,ashuffle);
    }

    free(bsh[0][0]);
    free(bsh[0]);
    free(bsh);
    free(bsh2[0][0]);
    free(bsh2[0]);
    free(bsh2);
    free(a1);
    free(a2);
    exit (0);
}

void coh1_init(int win, int m1, int m2, int l1, int l2)
/*< initialize >*/
{
	nw = win;
	n1 = m1;
	n2 = m2;
	lag1 = l1;
	lag2 = l2;
	u0 = sf_floatalloc3(2*nw+1, n1, n2);
	u1 = sf_floatalloc3(2*nw+1, n1, n2);
	v = sf_floatalloc2(2*(l1>l2? l1:l2)+1, n2);
	/* use two d memory for OPENMP */
}

int main(int argc, char* argv[]) 
{
    int wide1,wide2,wide3, wide, shift1, shift2, shift3, i, j, k, i1, n1, i2, n2, i3, n3, i4, n4;
    float ***data, ***signal, ***win;
    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_histint(in,"n2",&n2)) sf_error("No n2= in input");
    if (!sf_histint(in,"n3",&n3)) sf_error("No n3= in input");
    n4 = sf_leftsize(in,3);

    data = sf_floatalloc3(n1,n2,n3);
    signal = sf_floatalloc3(n1,n2,n3);

    if (!sf_getint("wide1",&wide1)) wide1=5;
    if (!sf_getint("wide2",&wide2)) wide2=5;
    if (!sf_getint("wide3",&wide3)) wide3=5;
    /* sliding window width */
    wide = wide1*wide2*wide3;

    win = sf_floatalloc3(wide1,wide2,wide3);

    for (i4=0; i4 < n4; i4++) {
	sf_floatread(data[0][0],n1*n2*n3,in);	
	for (i3=0; i3 < n3; i3++) { shift3 = SF_MAX (0, SF_MIN (n3-wide3, i3-wide3/2));
	    for (i2=0; i2 < n2; i2++) { shift2 = SF_MAX (0, SF_MIN (n2-wide2, i2-wide2/2));
		for (i1=0; i1 < n1; i1++) { shift1 = SF_MAX (0, SF_MIN (n1-wide1, i1-wide1/2));
		    for (i=0; i < wide; i++) {
			for (j=0; j < wide; j++) {
			    for (k=0; k < wide; k++) {
				win[i][j][k] = data[shift3+k][shift2+i][shift1+j];
			    }
			}
		    }
		    signal[i3][i2][i1] = sf_quantile(wide/2,wide,win[0][0]);
		}
	    }
	}
	
	sf_floatwrite(signal[0][0],n1*n2*n3,out);
    }	

    exit(0);
}

grid2 grid2_init (int n1, float o1, float d1 /* first axis */, 
		  int n2, float o2, float d2 /* second axis */,
		  float *slow2               /* data values [n1*n2] */, 
		  int order                  /* interpolation order */)
/*< Initialize grid object >*/
{
    int i, j;
    float v, vgrad[2];
    grid2 grd;

    grd = (grid2) sf_alloc(1,sizeof(*grd));

    grd->n1 = n1; grd->o1 = o1; grd->d1 = d1; 
    grd->n2 = n2; grd->o2 = o2; grd->d2 = d2;

    grd->pnt = sf_eno2_init (order != 0 ? order : 3, n1, n2);
    sf_eno2_set1 (grd->pnt, slow2);
    grd->vel = NULL;
    if (order <= 0) {
        grd->vel = sf_floatalloc3 (3, n1, n2);
        for (j = 0; j < n2; j++) {
            for (i = 0; i < n1; i++) {
                sf_eno2_apply (grd->pnt, i, j, 0., 0., &v, vgrad, BOTH);
                grd->vel[j][i][0] = slow2[j*n1 + i];
                grd->vel[j][i][1] = 0.5*vgrad[0]/d1;
                grd->vel[j][i][2] = 0.5*vgrad[1]/d2;
            }
        }
        sf_eno2_close (grd->pnt);
    }

    return grd;
}

void rgradient_init(char* type, int horder, int m1, int m2)
/*< initialize >*/
{
	float **p;
	int i;

	order = horder;
	nf = 2*order+1;
	c = sf_floatalloc2(nf, 2);
	p = lphpoly(order, order, type);
	for(i=0; i<2*nf; i++)	c[0][i] = p[0][i];
	free(p[0]);
	free(p);
	mode = type[0];

	n1 = m1;
	n2 = m2;

	b1 = sf_floatalloc2(n1, n2);
	b2 = sf_floatalloc2(n1, n2);

	g = sf_floatalloc3(n1*3, n2, nf);
	b = g[0];
	memset(g[0][0], 0, n1*n2*nf*sizeof(float));
#ifdef _OPENMP
    omp_init();
#endif
}

void* fcoh2_init(int n1, int n2, int pmin, int pmax, int nw)
/*< initialize >*/
{
    fcoh2 *p;

    p = sf_alloc(1,sizeof(fcoh2));
    p->n1 = n1;
    p->n2 = n2;
    p->min = pmin;
    p->max = pmax;
    p->ntw = nw;

    p->d1 = sf_floatalloc3(n1, n2, pmax-pmin+1);
    p->d2 = sf_floatalloc3(n1, n2, pmax-pmin+1);
    return p;
}

/*------------------------------------------------------------*/
slo3d slow3_init(cub3d   cub,
		 sf_fslice  slice_,   /* slowness slice */
		 int     nrmax,   /* maximum number of references */
		 float   dsmax,
		 float  twoway
    )
/*< initialize slowness >*/
{
    int imz, jj;
    int ompith=0;

    /*------------------------------------------------------------*/
    slo3d slo;
    slo = (slo3d) sf_alloc(1,sizeof(*slo));

    slo->slice=slice_;
    slo->nrmax=nrmax;
    slo->dsmax=dsmax;

    if(twoway) { slo->twoway = 2;
    } else {     slo->twoway = 1;
    }

    slo->ss = sf_floatalloc3(cub->alx.n,cub->aly.n,cub->ompnth);  /* slowness */
    slo->so = sf_floatalloc3(cub->alx.n,cub->aly.n,cub->ompnth);  /* slowness */
    slo->sm = sf_floatalloc2(slo->nrmax,cub->amz.n);  /* ref slowness squared */
    slo->nr = sf_intalloc              (cub->amz.n);  /* nr of ref slownesses */
    
    for (imz=0; imz<cub->amz.n; imz++) {
	sf_fslice_get(slo->slice,imz,slo->ss[0][0]);
	slow3_twoway(cub,slo,slo->ss,ompith);

	slo->nr[imz] = slow3(slo->nrmax,
			     slo->dsmax,
			     cub->alx.n*cub->aly.n,
			     slo->ss[0][0],
			     slo->sm[imz]);
    }
    for (imz=0; imz<cub->amz.n-1; imz++) {
	for (jj=0; jj<slo->nr[imz]; jj++) {
	    slo->sm[imz][jj] = 0.5*(slo->sm[imz][jj]+slo->sm[imz+1][jj]);
	}
    }

    return slo;
}

void resample_vel(const modeling_t *olds, const modeling_t *news, const vel_t *oldv, vel_t *newv) {
  free(**newv->dat); free(*newv->dat); free(newv->dat);
  free(**newv->vgamma); free(*newv->vgamma); free(newv->vgamma);

  newv->dat = sf_floatalloc3(news->n1, news->n2, news->n3);
  newv->vgamma = sf_floatalloc3(news->n1, news->n2, news->n3);

  interpfield(olds, news, oldv->dat, newv->dat, true);

  for (int i3 = 0; i3 < news->n3; i3++) {
    for (int i2 = 0; i2 < news->n2; i2++) {
      for (int i1 = 0; i1 < news->n1; i1++) {
        newv->vgamma[i3][i2][i1] = -powf(newv->dat[i3][i2][i1],2*gs_gamma-1) * powf(gs_w0, 2 * gs_gamma) * sin(SF_PI * gs_gamma) / news->dt;
      }
    }
  }
}

void resample_p(const modeling_t *olds, const modeling_t *news, float ****p)
{
  float ***newp = sf_floatalloc3(news->n1, news->n2, news->n3);
  interpfield(olds, news, *p, newp, false);

  free((*p)[0][0]); free((*p)[0]); free(*p);
  *p = newp;
}

int main(int argc, char* argv[])
{
    int n1, n2,o1,d1,d2,i1,i2,ix,iz; 

    float  ***dat,**sum;
    float sum1; 
    sf_file in, out; /* Input and output files */

    /* Initialize RSF */
    sf_init(argc,argv);
    /* standard input */
    in  = sf_input("in");
    /* standard output */
    out = sf_output("out");
    
    /* parameters from input file*/
    if (SF_FLOAT != sf_gettype(in)) sf_error("Need float input");
    if (!sf_histint(in,"o1",&o1)) sf_error("No o1= in input");
    if (!sf_histint(in,"d1",&d1)) sf_error("No d1= in input");
    if (!sf_histint(in,"d2",&d2)) sf_error("No d2= in 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");
    if (!sf_histint(in,"n3",&n3)) sf_error("No n3= in input");
    fprintf(stderr,"values: %d %d\n",n1,n2,n3);
   sum = sf_floatalloc2(n3,n3);
     dat = sf_floatalloc3 (n1,n2,n3);

    fprintf(stderr,"values: %d %d\n",n1,n2);
  dat = sf_floatalloc2 (n1,n2);

   for (ix=0; ix<n2; ++ix) {
      for (iz=0; iz<n1; ++iz) {

           dat[ix][iz]=0.0;
      }
     }
sum1=0;

    fprintf(stderr,"values1: %f\n",sum1);
    sf_floatread(dat[0],n1*n2,in); 

//sum1=0;

    fprintf(stderr,"values1: %f\n",sum1);
      for (i2=o1; i2<n2; ++i2){
            for (i1=o1; i1<n1; ++i1){
      sum1 += dat[i2][i1];
        }
        }
    fprintf(stderr,"values: %f\n",sum1);
   sum[0][0]=sum1; 
    sf_floatwrite(sum[0],1*1,out);
      

  exit(0);
}

vel_t *clone_vel(float ***vel, int nz, int nx, int ny,
    float oz, float ox, float oy,
    float dz, float dx, float dy,
    float w0, float qfact)
{
  vel_t *v = malloc(sizeof *v);
  v->n1 = nz; v->n2 = nx; v->n3 = ny;
  v->o1 = oz; v->o2 = ox; v->o3 = oy;
  v->d1 = dz; v->d2 = dx; v->d3 = dy;
  v->dat = sf_floatalloc3(nz, nx, ny);
  memcpy(v->dat[0][0], vel[0][0], nz*nx*ny*sizeof(float));
  v->vgamma = sf_floatalloc3(1,1,1);

  gs_w0 = w0; gs_qfact= qfact;
  gs_gamma = 1.0/SF_PI*atan(2*SF_PI/qfact);

  sf_warning("gs_gamma: %f", gs_gamma);
  return v;
}

/*------------------------------------------------------------*/
wexop3d wex_init(wexcub3d cub)
/*< initialize >*/
{
    wexop3d weop;
    weop = (wexop3d) sf_alloc(1,sizeof(*weop));

    weop->wws = sf_complexalloc3(cub->amx.n,cub->amy.n,cub->ompnth);
    weop->wwr = sf_complexalloc3(cub->amx.n,cub->amy.n,cub->ompnth);
    weop->ws  = sf_complexalloc4(cub->amx.n,cub->amy.n,cub->az.n,cub->ompnth);
    weop->wr  = sf_complexalloc4(cub->amx.n,cub->amy.n,cub->az.n,cub->ompnth);

    weop->rr  = sf_floatalloc3(cub->amx.n,cub->amy.n,cub->az.n);

    return weop;
}

/*------------------------------------------------------------*/
void fdbell3d_init(int n)
/*< init bell taper >*/
{
    int   i1,i2,i3;
    float s;

    nbell = n;
    s = 0.5*nbell;

    bell3d=sf_floatalloc3(2*nbell+1,2*nbell+1,2*nbell+1);

    for        (i3=-nbell;i3<=nbell;i3++) {
	for    (i2=-nbell;i2<=nbell;i2++) {
	    for(i1=-nbell;i1<=nbell;i1++) {
		bell3d[nbell+i3][nbell+i2][nbell+i1] = exp(-(i1*i1+i2*i2+i3*i3)/s);
	    }
	}    
    }
}

void weimc_init(bool verb_,
		axa amx_        /* i-line (data) */,
		axa amy_        /* x-line (data) */,
		axa amz_        /* depth */,
		axa aw_         /* frequency */
)
/*< initialize WE IC >*/
{
    verb=verb_;

    amz = amz_;
    amx = amx_;
    amy = amy_;
    aw  = aw_;

    /* allocate wavefield storage */
    us = sf_complexalloc3(amx.n,amy.n,amz.n);
    ur = sf_complexalloc3(amx.n,amy.n,amz.n);
    qq = sf_floatalloc3  (amx.n,amy.n,amz.n);
}