C++ NINT函数代码示例

void rc1_fft(REAL *data, complex *cdata, int n, int sign)
{
	int    j;
	double *datft;

	if (NINT(pow(2.0, (double)NINT(log((double)n)/log(2.0)))) != n) {
		if (npfar(n) == n) pfarc(sign,n,data,cdata);
		else rcdft(data,cdata,n,sign);
	}
	else {
		datft = (double *)malloc(n*sizeof(double));
		if (datft == NULL) fprintf(stderr,"rc1_fft: memory allocation error\n");
	
		for (j = 0; j < n; j++) datft[j] = (double)data[j];
		realfft(n, datft);
		cdata[0].i = 0.0;
		for (j = 0; j < n/2; j++) {
			cdata[j].r = (REAL)datft[j]; 
			cdata[j+1].i = sign*(REAL)datft[n-j-1]; 
		}
		cdata[n/2].r = datft[n/2];
		cdata[n/2].i = 0.0; 
	
		free(datft);
	}

	return;
}

void cr1_fft(complex *cdata, REAL *data, int n, int sign)
{
	int    j;
	double *datft;

	if (NINT(pow(2.0, (double)NINT(log((double)n)/log(2.0)))) != n) {
		if (npfar(n) == n) pfacr(sign,n,cdata,data);
		else crdft(cdata,data,n,sign);
	}
	else {
		datft = (double *)malloc(n*sizeof(double));
		if (datft == NULL) fprintf(stderr,"cr1_fft: memory allocation error\n");

		for (j = 0; j < n/2; j++) {
			datft[j] = (double)cdata[j].r;
			datft[n-1-j] = (double)cdata[j+1].i;
		}
		datft[n/2] = (double)cdata[n/2].r;

		realifft(n, datft);
	
		if (sign == -1) {
			for (j = 0; j < n; j++) data[j] = (REAL)datft[j];
		}
		else if (sign == 1) {
			for (j = 1; j < n; j++) data[j] = (REAL)datft[n-j];
			data[0] = (REAL)datft[0];
		}
	
		free(datft);
	}
	
	return;
}

int FromLT (int rsize, double *ltm, double *ltv, int *bin, int *corner) {

/* arguments:
int rsize        i: reference pixel size
double ltm[2]    i: diagonal elements of MWCS matrix
double ltv[2]    i: MWCS linear transformation vector
int bin[2]       o: pixel size in X and Y
int corner[2]    o: corner of subarray in X and Y
*/

	extern int status;
	double dbinx, dbiny, dxcorner, dycorner;

	dbinx = (double)rsize / ltm[0];
	dbiny = (double)rsize / ltm[1];

	dxcorner = (dbinx - rsize) - dbinx * ltv[0];
	dycorner = (dbiny - rsize) - dbiny * ltv[1];

	/* Round off to the nearest integer. */
	corner[0] = NINT (dxcorner);
	corner[1] = NINT (dycorner);
	bin[0] = NINT (dbinx);
	bin[1] = NINT (dbiny);

	return (status);
}

asynStatus MCB4BAxis::sendAccelAndVelocity(double acceleration, double velocity) 
{
  asynStatus status;
  int ival;
  // static const char *functionName = "MCB4B::sendAccelAndVelocity";

  // Send the velocity
  ival = NINT(fabs(115200./velocity));
  if (ival < 2) ival=2;
  if (ival > 255) ival = 255;
  sprintf(pC_->outString_, "#%02dV=%d", axisNo_, ival);
  status = pC_->writeReadController();

  // Send the acceleration
  // acceleration is in steps/sec/sec
  // MCB is programmed with Ramp Index (R) where:
  //   dval (steps/sec/sec) = 720,000/(256-R) */
  //   or R=256-(720,000/dval) */
  ival = NINT(256-(720000./acceleration));
  if (ival < 1) ival=1;
  if (ival > 255) ival=255;
  sprintf(pC_->outString_, "#%02dR=%d", axisNo_, ival);
  status = pC_->writeReadController();
  return status;
}

void cc1_fft(complex *cdata, int n, int sign)
{
	int    j;
	double  *real, *imag;

	if (NINT(pow(2.0, (double)NINT(log((double)n)/log(2.0)))) != n) {
		if (npfa(n) == n) pfacc(sign, n, cdata);
		else ccdft(cdata,n,sign);
	}
	else {
		real = (double *)malloc(n*sizeof(double));
		if (real == NULL) fprintf(stderr,"cc1_fft: memory allocation error\n");
		imag = (double *)malloc(n*sizeof(double));
		if (imag == NULL) fprintf(stderr,"cc1_fft: memory allocation error\n");
	
		for (j = 0; j < n; j++) {
			real[j] = (double)cdata[j].r;
			imag[j] = (double)cdata[j].i;
		}

		if (sign < 0) fft(n, real, imag);
		else ifft(n, real, imag);

		for (j = 0; j < n; j++) {
			cdata[j].r = (REAL)real[j];
			cdata[j].i = (REAL)imag[j];
		}

		free(real);
		free(imag);
	}

	return;
}

void smooth (float **in, float **out, long naxes[], float sigma, float radius)
{
    int i, j, k, l, xmin, ymin, xmax, ymax;
    float sumwt, sumflux, wt, rad;

    for (j=1; j <= naxes[2]; j++) {
	for (i=1; i <= naxes[1]; i++) {
	    xmin = IMAX (1, i - NINT(radius));
	    xmax = IMIN (i+NINT(radius), naxes[1]);
	    ymin = IMAX (1, j - NINT(radius));
	    ymax = IMIN (j+NINT(radius), naxes[2]);
	    sumwt = 0.;
	    sumflux = 0.;
	    for (l=ymin; l <= ymax; l++) {
		for (k=xmin; k <= xmax; k++) {
		    rad = sqrt((k-i)*(k-i)+(l-j)*(l-j));
		    if (rad <= radius) {
		        wt = exp (-rad*rad/2/sigma/sigma);
			sumwt += wt;
			sumflux += in[l][k] * wt;
		    };
		};	    
	    };
	    out[j][i] = sumflux / sumwt;
	};
    };
}

void makefilter(float *f, int nfft, int nfreq, float dt, float *filter)
{
        float onfft = 1.0 / nfft;
        float df = onfft / dt;
        int nfreqm1 = nfreq - 1;
        int if1 = NINT(f[0]/df);
        int if2 = NINT(f[1]/df);
        int if3 = MIN(NINT(f[2]/df), nfreqm1);
        int if4 = MIN(NINT(f[3]/df), nfreqm1);


        { register int i;
	  register float c = PIBY2 / (if2 - if1 + 2);
	  for (i = if1; i <= if2; ++i) {
		register float s = sin(c*(i - if1 + 1));
		filter[i] = s * s * onfft;
          }
	 }

        { register int i;
	  register float c = PIBY2 / (if4 - if3 + 2);
	  for (i = if3; i <= if4; ++i) {
		register float s = sin(c*(if4 - i + 1));
		filter[i] = s * s * onfft;
	  }
        }

        { register int i;
          for (i = if2 + 1; i < if3;   ++i)  filter[i] = onfft; 
          for (i = 0;       i < if1;   ++i)  filter[i] = 0.0; 
          for (i = if4 + 1; i < nfreq; ++i)  filter[i] = 0.0; 
        }
}

/*==========================================================================*/
void lnk_cell_dis(int ia, int ib, int ic, 
		  int ncell_a, int ncell_b, int ncell_c, 
		  double *rcell, double hmat_lnk[], int iperd)

/*==========================================================================*/
{/*Begin Routine */
/*==========================================================================*/
/*                 Local variables                                         */
  double acelli, bcelli, ccelli, dx, dy, dz, dsx, dsy, dsz;
  double dsxt, dsyt, dszt;

/*==========================================================================*/
/*  ia,ib,ic are the shifts along a,b,c                                     */
/*  ncell_a,ncell_b,ncell_c are the number of divisions along a,b,c         */
/*  hmat matrix of cell parameters                                          */
/*  rcell is return with cartesian distance associated with shifts          */
/* scalar temp                                                              */
/* get inverse number of divisions along cell axis                          */
/*==========================================================================*/

  acelli = 1. / (double) ncell_a;
  bcelli = 1. / (double) ncell_b;
  ccelli = 1. / (double) ncell_c;

  /* convert shifts to a displacement along fractional coord. (the s) */

  dsx = (double) ia * acelli;
  dsy = (double) ib * bcelli;
  dsz = (double) ic * ccelli;

  /* periodic image */
/*  if (iperd >= 1) {  dsx -= NINT(dsx);}*/
/*  if (iperd >= 2) {  dsy -= NINT(dsy);}*/
/*  if (iperd >= 3) {  dsz -= NINT(dsz);}*/
    dsx -= NINT(dsx);
    dsy -= NINT(dsy);
    dsz -= NINT(dsz);

  /* subtract out one cell to account that cell that share edges can have */
  /* particles separteded by one less cell distances */

/*  dsxt = copysign(acelli,dsx); */
/*  dsyt = copysign(bcelli,dsy); */
/*  dszt = copysign(ccelli,dsz); */
  if(dsx>=0){dsxt=acelli;}else{dsxt=-acelli;}
  if(dsy>=0){dsyt=bcelli;}else{dsyt=-bcelli;}
  if(dsz>=0){dszt=ccelli;}else{dszt=-ccelli;}
  if (dsx != 0.) {dsx -= dsxt;}
  if (dsy != 0.) {dsy -= dsyt;}
  if (dsz != 0.) {dsz -= dszt;}

  /* convert to cartesian displacements */
  dx = dsx * hmat_lnk[1] + dsy * hmat_lnk[4] + dsz * hmat_lnk[7];
  dy = dsx * hmat_lnk[2] + dsy * hmat_lnk[5] + dsz * hmat_lnk[8];
  dz = dsx * hmat_lnk[3] + dsy * hmat_lnk[6] + dsz * hmat_lnk[9];

  /* get distance */
  *rcell = sqrt(dx * dx + dy * dy + dz * dz);
} /* lnk_cell_dis */

int
main(int argc, char **argv)
{
	int i, j, n;
	float *f,  max, step, rstep;
	float fhist[1024], dev, rdev, ent, error; 
	int hist[1024];
	 
	initargs(argc, argv);
	requestdoc(1);

	MUSTGETPARINT("n",&n);

	f = alloc1float(n);

	fread(f,sizeof(float),n,stdin);
	
	for(i=0;i<1024;i++) hist[i] = 0;

	for(i=0,rdev=0.;i<n;i++)
	   rdev += f[i]*f[i]; 
	rdev = rdev/n;
	rdev = sqrt(rdev);

	if(!getparfloat("dev",&dev)) dev = rdev;

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

	step = dev*3.464*.01;

	rstep = 1./step;

	error = 0.;
	for(i=0;i<n;i++){
	    max = f[i]*rstep;
	    error += (NINT(max)*step - f[i])*(NINT(max)*step - f[i]);
	    hist[NINT(max)+512] ++;
	}

	error = error/n;
	error = sqrt(error);
	error = error/rdev;

	ent = 0.;
	for(j=0;j<1024;j++){
		fhist[j] = ((float) hist[j])/((float) n);
		if(hist[j])
		    ent += fhist[j]*log(fhist[j])/log(2.);
	}
	ent = -ent;

fprintf(stderr,"entropy of the signal is=%f, average error=%f\n",ent, error);

	fwrite(fhist,sizeof(float),1024,stdout);

	return EXIT_SUCCESS;
}

void remove_fb(float *yp,float *ym,int n,short *scaler,short *shft)
/* Find Scale and timeshift 

	Yp = data trace
	Ym = wavelet
	
	F=min(Yp(x) - Ym(x)*a)^2 is a function to minimize for scaler
	
	find shift first with xcorrelation then
	
	solve for a - scaler 
*/



{

#define RAMPR 5

	void find_p(float *ym,float *yp,float *a,int *b,int n);
	
	int it,ir;
	
	float a=1.0;
	int b=0;
	int ramps;
	float *a_ramp;
	
	ramps=NINT(n-n/RAMPR);
	
	find_p(ym,yp,&a,&b,n);
	
	a_ramp = ealloc1float(n);
	
	for(it=0;it<ramps;it++)
		a_ramp[it]=1.0;
	
	for(it=ramps,ir=0;it<n;it++,ir++) {
		a_ramp[it]=1.0-(float)ir/(float)(n-ramps-1);
	}
	
	
	
	if (b<0) {
		for(it=0;it<n+b;it++) 
			yp[it-b] -=ym[it]*a*a_ramp[it-b];
	} else {
		for(it=0;it<n-b;it++) 
			yp[it] -=ym[it-b]*a*a_ramp[it+b];
	}
	
	*scaler = NINT((1.0-a)*100.0);
	*shft = b;
	
	free1float(a_ramp);
}

/*------------------------------------------------------------*/
lint2d lint2d_make(	int na, 
					pt2d *aa, 
					fdm2d fdm)
		   
/*< init 2D linear interpolation >*/
{
    lint2d ca;
    int    ia;
    float f1,f2;
    
    ca = (lint2d) sf_alloc(1,sizeof(*ca));

    ca->n = na;

    ca->w00 = sf_floatalloc(na);
    ca->w01 = sf_floatalloc(na);
    ca->w10 = sf_floatalloc(na);
    ca->w11 = sf_floatalloc(na);

    ca->jz  = sf_intalloc(na);
    ca->jx  = sf_intalloc(na);

    for (ia=0; ia<na; ia++) {
	
		if(aa[ia].z >= fdm->ozpad && 
			aa[ia].z <  fdm->ozpad + (fdm->nzpad-1)*fdm->dz &&
			aa[ia].x >= fdm->oxpad && 
			aa[ia].x <  fdm->oxpad + (fdm->nxpad-1)*fdm->dx   ) {
				
			ca->jz[ia] = NINT((aa[ia].z-fdm->ozpad)/fdm->dz);
			ca->jx[ia] = NINT((aa[ia].x-fdm->oxpad)/fdm->dx);
				
			f1 = (aa[ia].z-fdm->ozpad)/fdm->dz - ca->jz[ia];
			f2 = (aa[ia].x-fdm->oxpad)/fdm->dx - ca->jx[ia];
		}
		else {
			ca->jz[ia] = 0; 
			ca->jx[ia] = 0;
			
			f1 = 1; 
			f2 = 0;
		}
		
		ca->w00[ia] = (1-f1)*(1-f2);
		ca->w01[ia] = (  f1)*(1-f2);
		ca->w10[ia] = (1-f1)*(  f2);
		ca->w11[ia] = (  f1)*(  f2);
		
    }
	
    return ca;
}

static void FirstLast (double *ltm, double *ltv, int *snpix, int *npix,
		int *rbin, int *first, int *last, int *sfirst) {

/* arguments:
double ltm[2], ltv[2]   i: linear transformation from scratch to image coords
				(adjusted for use with zero indexed coords)
int snpix[2]            i: size of scratch image
int npix[2]             i: size of actual image
int rbin[2]             o: number of pixels in scratch for one image pixel
int first[2], last[2]   o: corners of overlap region, in image coords
int sfirst[2]           o: lower left corner of overlap, in scratch array
*/

	double scr;		/* pixel coordinate in scratch array */
	int i, k;

	rbin[0] = NINT (1. / ltm[0]);
	rbin[1] = NINT (1. / ltm[1]);

	for (k = 0;  k < 2;  k++) {	/* axis number */
	    /* Search for the first pixel in the image array that maps to a
		point that is completely within the scratch array.  The left
		(lower) edge of pixel i is (i - 0.5).  Map (i - 0.5) to the
		scratch array, and if that point is within the scratch array,
		i is the first fully illuminated pixel.
	    */
	    for (i = 0;  i < npix[k];  i++) {
		scr = (i - 0.5 - ltv[k]) / ltm[k];
		/* -0.5 is left (lower) edge of first pixel in scratch */
		if (scr+TOLERANCE >= -0.5) {
		    first[k] = i;
		    sfirst[k] = NINT (scr+0.5);
		    break;
		}
	    }
	    /* Now look for the last fully illuminated pixel, using the
		right (upper) edge of the image pixels.
	    */
	    for (i = npix[k]-1;  i > 0;  i--) {
		scr = (i + 0.5 - ltv[k]) / ltm[k];
		/* compare scr with right (upper) edge of last pixel in the
		   scratch array
		*/
		if (scr-TOLERANCE <= snpix[k]-1. + 0.5) {
		    last[k] = i;
		    break;
		}
	    }
	}
}

void crm_fft(complex *cdata, REAL *data, int n1, int n2, int ldc, int ldr, int sign)
{
	int    j, i;
	double *datft;

	if (NINT(pow(2.0, (double)NINT(log((double)n1)/log(2.0)))) != n1) {
		if (npfar(n1) == n1) {
			if (ldr == n1 && ldc == n2) {
				pfa2cr(sign, 1, n1, n2, cdata, data);
			}
			else {
				for (i = 0; i < n2; i++) {
					pfacr(sign, n1, &cdata[i*ldc], &data[i*ldr]);
				}
			}
		}
		else {
			for (i = 0; i < n2; i++) {
				crdft(&cdata[i*ldc], &data[i*ldr], n1, sign);
			}
		}
	}
	else {
		datft = (double *)malloc(n1*sizeof(double));
		if (datft == NULL) fprintf(stderr,"crm_fft: memory allocation error\n");
	
		for (i = 0; i < n2; i++) {
			for (j = 0; j < n1/2; j++) {
				datft[j] = (double)cdata[i*ldc+j].r;
				datft[n1-1-j] = (double)cdata[i*ldc+j+1].i;
			}
			datft[n1/2] = (double)cdata[i*ldc+n1/2].r;
	
			realifft(n1, datft);
	
			if (sign == -1) {
				for (j = 0; j < n1; j++) data[i*ldr+j] = (REAL)datft[j];
			}
			else if (sign == 1) {
				for (j = 1; j < n1; j++) data[i*ldr+j] = (REAL)datft[n1-j];
				data[i*ldr] = (REAL)datft[0];
			}
		}
	
		free(datft);
	}

	return;
}

asynStatus MCB4BAxis::move(double position, int relative, double minVelocity, double maxVelocity, double acceleration)
{
  asynStatus status;
  // static const char *functionName = "MCB4BAxis::move";

  status = sendAccelAndVelocity(acceleration, maxVelocity);
  
  if (relative) {
    sprintf(pC_->outString_, "#%02dI%+d", axisNo_, NINT(position));
  } else {
    sprintf(pC_->outString_, "#%02dG%+d", axisNo_, NINT(position));
  }
  status = pC_->writeReadController();
  return status;
}

static void adjustAxesValues (int n1, float d1, float f1,
	int n2, float d2, float f2,
	float *x1beg, float *x1end, float *x2beg, float *x2end)
/*****************************************************************************
Adjust axes values to nearest samples; ensure at least 2 by 2 samples.
*****************************************************************************/
{
	int ix1beg,ix1end,ix2beg,ix2end;
	
	ix1beg = NINT((*x1beg-f1)/d1);
	ix1end = NINT((*x1end-f1)/d1);
	ix1beg = MAX(0,MIN(n1-1,ix1beg));
	ix1end = MAX(0,MIN(n1-1,ix1end));
	if (ix1beg==ix1end) {
		if (*x1beg<*x1end) {
			if (ix1beg>0)
				ix1beg--;
			else
				ix1end++;
		} else {
			if (ix1beg<n1-1)
				ix1beg++;
			else
				ix1end--;
		}
	}
	*x1beg = f1+ix1beg*d1;
	*x1end = f1+ix1end*d1;
	ix2beg = NINT((*x2beg-f2)/d2);
	ix2end = NINT((*x2end-f2)/d2);
	ix2beg = MAX(0,MIN(n2-1,ix2beg));
	ix2end = MAX(0,MIN(n2-1,ix2end));
	if (ix2beg==ix2end) {
		if (*x2beg<*x2end) {
			if (ix2beg>0)
				ix2beg--;
			else
				ix2end++;
		} else {
			if (ix2beg<n2-1)
				ix2beg++;
			else
				ix2end--;
		}
	}	
	*x2beg = f2+ix2beg*d2;
	*x2end = f2+ix2end*d2;
}