C++ DoubleArray类代码示例

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor::computeNoise(const DoubleArray &data, int i1, int i2, double *offset, double *amplitude) {
	// compute offset and rms within the time window
	if(i1<0) i1=0;
	if(i2<0) return false;
	if(i2>(int)data.size()) i2=(int)data.size();

	// If noise window is zero return an amplitude and offset of zero as well.
	if ( i2-i1 == 0 ) {
		*amplitude = 0;
		*offset = 0;
		return true;
	}

	DoubleArrayPtr d = static_cast<DoubleArray*>(data.slice(i1, i2));

	double ofs, amp;

	// compute pre-arrival offset
	ofs = d->median();
	// compute rms after removing offset
	amp = 2 * d->rms(ofs);

	if ( offset ) *offset = ofs;
	if ( amplitude ) *amplitude = amp;

	return true;
}

int main(int argc, char* argv[])
{

	ifstream in("../dict/words.lex");                 // input

	if(!in) {
		cout << "Cannot open input file.\n";
		return 1;
	}

	char ** buf = new char*[655350];
	int i=0;
	while(in) {
		char * str = new char[64];
		in.getline(str, 64);                // delim defaults to '\n'
		//if(in) 
		//	cout << str << endl;
		buf[i] = str;
		i++;
		if (i % 100 == 0)
			cout << i << str << endl;
	}

	in.close();

	DoubleArray* da = new DoubleArray();
	da->build(i, (char**)buf);

	da->save("../dict/dartswords.lex");

	return 0;
}

/*
This has a funny way of dealing with edges. In case of a Gaussian blur, divides by the sum of all of the numbers in the mask by definition. Might prove harsh on edges. In all other cases, divides by 8 no matter what, so that can equalize at edges. Also, when one of the pixels needed in the mask calculation is off the bounds of the image, the center pixel value is substituted, so that there is less difference. This doesnt pick up as many stray lines on edges, but also proves to often not pick up enough lines. 
*/
DoubleArray Mask::applyMask(DoubleArray array){
	ICoord size = array.size();
	DoubleArray newarray = DoubleArray(size);
	width = (int)maskArray.size()(0)/2; /* in case the maskArray was edited */
	for (DoubleArray::iterator i = array.begin();i !=array.end(); ++i){
		ICoord curr = i.coord();
		double num = 0;
		double c = 0;
		int d = 1;
		for (int i = -width; i <= width; ++i){ // scans through the mask
			for (int j = -width; j <= width; ++j){
				ICoord a = ICoord(i,j);
				if (pixelInBounds(curr + a, size)){
					num = num + array[curr+a]*maskArray[ICoord(i + width, j + width)];
					c = c + maskArray[ICoord(i + width, j + width)];
				}
				else{
					num = num + array[curr]*maskArray[ICoord(i + width, j + width)]; /* if pixel needed for mask is out of the image, then substitute the center pixel */
					d = d + 1; /* count number of pixels out of bounds */
				}
			}
		}
		if (type == GAUSSIAN_MASK)
			newarray[curr] = num/159;
		else if (type == SMALL_GAUSSIAN_MASK)
			newarray[curr] = num/99;
		else
			newarray[curr] = num/8;
	}
	return newarray;
}

/*!
	Updates the ImageView and the ParameterSelector associated
	with the given panel.
*/
void VideoParserWindow::UpdatePanel(VideoControlPanel* pPanel, int sel, 
		DoubleArray params, bool checkViewSync)
{
	// Get the potentially new parameter range from the component
	DoubleArray minVals, maxVals, steps;

	m_videoProcessor.GetParameterInfo(sel, &minVals, &maxVals, &steps);

	// Ensure that the current parameters are within the "new" ranges
	if (params.size() == minVals.size())
	{
		for (unsigned i = 0; i < params.size(); i++)
		{
			if (params[i] > maxVals[i])
				params[i] = maxVals[i];

			if (params[i] < minVals[i])
				params[i] = minVals[i];
		}
	}
	else
	{
		params = minVals;
	}

	// Update the selected view with the new image
	UpdateImageView(pPanel, sel, params, checkViewSync);

	// Update the displayed parameter value and range
	pPanel->GetParamSelector()->Update(params, minVals, maxVals, steps);

	pPanel->GetCommandSelector()->UpdateCommands(
		m_videoProcessor.GetUserCommands(sel));
}

	int Ardb::ZAdd(const DBID& db, const Slice& key, double score,
			const Slice& value)
	{
		DoubleArray scores;
		SliceArray svs;
		scores.push_back(score);
		svs.push_back(value);
		return ZAdd(db, key, scores, svs);
	}

//--------------------------------------------------------------------------------------------------
/// 
//--------------------------------------------------------------------------------------------------
TEST(ArrayTest, CopyConvertedData)
{
    // Double array to float
    {
        DoubleArray ad;
        ad.resize(4);
        ad[0] = 0.0;
        ad[1] = 1.0;
        ad[2] = 2.0;
        ad[3] = 3.0;

        // Copy full array
        FloatArray af;
        af.resize(4);
        af.copyConvertedData(ad, 4, 0, 0);
        EXPECT_FLOAT_EQ(0.0f, af[0]);
        EXPECT_FLOAT_EQ(1.0f, af[1]);
        EXPECT_FLOAT_EQ(2.0f, af[2]);
        EXPECT_FLOAT_EQ(3.0f, af[3]);

        // Copy partial array to float array
        af.resize(2);
        af.setAll(0);
        af.copyConvertedData(ad, 2, 0, 1);

        EXPECT_FLOAT_EQ(1.0f, af[0]);
        EXPECT_FLOAT_EQ(2.0f, af[1]);
    }

    // Vec3d to Vec3f and Vec3i
    {
        Vec3dArray ad;
        ad.resize(2);
        ad[0].set(1.1, 2.5, 3.9);
        ad[1].set(11.1, 12.5, 13.9);

        Vec3fArray af;
        af.resize(2);
        af.copyConvertedData(ad, 2, 0, 0);
        EXPECT_FLOAT_EQ(1.1f,  af[0].x());
        EXPECT_FLOAT_EQ(2.5f,  af[0].y());
        EXPECT_FLOAT_EQ(3.9f,  af[0].z());
        EXPECT_FLOAT_EQ(11.1f, af[1].x());
        EXPECT_FLOAT_EQ(12.5f, af[1].y());
        EXPECT_FLOAT_EQ(13.9f, af[1].z());

        Array<Vec3i> ai;
        ai.resize(2);
        ai.copyConvertedData(ad, 2, 0, 0);
        EXPECT_EQ(1,  ai[0].x());
        EXPECT_EQ(2,  ai[0].y());
        EXPECT_EQ(3,  ai[0].z());
        EXPECT_EQ(11, ai[1].x());
        EXPECT_EQ(12, ai[1].y());
        EXPECT_EQ(13, ai[1].z());
    }
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor_msbb::computeAmplitude(const DoubleArray &data,
                                               size_t i1, size_t i2,
                                               size_t si1, size_t si2, double offset,
                                               AmplitudeIndex *dt,
                                               AmplitudeValue *amplitude,
                                               double *period, double *snr) {
	/*
	* Low-level signal amplitude computation. This is magnitude specific.
	*
	* Input:
	*      f           double array of length n
	*      i1,i2       indices defining the measurement window,
	*                  0 <= i1 < i2 <= n
	*      offset      this is subtracted from the samples in f before
	*                  computation
	*
	* Output:
	*      dt          Point at which the measurement was mad/completed. May
	*                  be the peak time or end of integration.
	*      amplitude   amplitude. This may be a peak amplitude as well as a
	*                  sum or integral.
	*      period      dominant period of the signal. Optional. If the period
	*                  is not computed, set it to -1.
	*/

	size_t imax = find_absmax(data.size(), (const double*)data.data(), si1, si2, offset);
	double amax = fabs(data[imax] - offset);
	double pmax = -1;
	double pstd =  0; // standard error of period
	if ( !measure_period(data.size(), static_cast<const double*>(data.data()), imax, offset, &pmax, &pstd) )
		pmax = -1;

	if ( amax < *_noiseAmplitude * _config.snrMin ) {
		amplitude->value = amax / *_noiseAmplitude;
		setStatus(LowSNR, amplitude->value);
		return false;
	}

	dt->index = imax;
	*period = pmax;
	amplitude->value = amax;

	if ( _usedComponent <= SecondHorizontal ) {
		if ( _streamConfig[_usedComponent].gain != 0.0 )
			amplitude->value /= _streamConfig[_usedComponent].gain;
		else {
			setStatus(MissingGain, 0.0);
			return false;
		}
	}
	else
		return false;

	return true;
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor::deconvolveData(Response *resp, DoubleArray &data,
                                        int numberOfIntegrations) {
	// Remove linear trend
	double m,n;
	Math::Statistics::computeLinearTrend(data.size(), data.typedData(), m, n);
	Math::Statistics::detrend(data.size(), data.typedData(), m, n);

	return resp->deconvolveFFT(data, _stream.fsamp, _config.respTaper,
	                           _config.respMinFreq, _config.respMaxFreq,
	                           numberOfIntegrations);
}

void doClassifying(DoubleArray & gray) {
    ICoord size = gray.size();
    DoubleArray newgray = DoubleArray(size);
    ClassifyClass classif = ClassifyClass();
    newgray = classif.classifyRegions(gray);
    gray = newgray;

}

DoubleArray MASKS::applyMasks(const DoubleArray& image) {
//This applies the mask that is stored onto an array of doubles representing
//the image.
  DoubleArray newImage(image.size(),0.0);
  double tempSum=0;
  int startI,endI,startK,endK;
  int M=maskArrays.width(), N=maskArrays.height();

  for(int x=0;x<image.width();x++) {
    for(int y=0;y<image.height();y++) {
      startI=-M/2; endI=M/2;
      startK=-N/2; endK=N/2;
      if(x-M/2<0)
	startI=-x;
      if(y-N/2<0)
	startK=-y;
      if(x+M/2>=image.width())
	endI=image.width()-x-1;
      if(y+N/2>=image.height())
	endK=image.height()-y-1;
      for(int i=startI;i<=endI;i++) {
	for(int k=startK;k<=endK;k++) {
	  tempSum+=image[ICoord(x+i,y+k)]*maskArrays[ICoord(M/2+i,N/2+k)];
	  //If M and N aren't odd, will get segmentation fault
	}
      }
      newImage[ICoord(x,y)]=tempSum;
      //    if(tempSum<0)
      tempSum=0;
    }
  }
  return newImage;
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor_Mjma::computeAmplitude(
		const DoubleArray &data,
		size_t i1, size_t i2,
		size_t si1, size_t si2,
		double offset,AmplitudeIndex *dt,
		AmplitudeValue *amplitude,
		double *period, double *snr)
{
	double amax;

	size_t imax = find_absmax(data.size(), data.typedData(), si1, si2, offset);
	amax = fabs(data[imax] - offset);
	dt->index = imax;

	if ( *_noiseAmplitude == 0. )
		*snr = 1000000.0;
	else
		*snr = amax / *_noiseAmplitude;

	if ( *snr < _config.snrMin ) {
		setStatus(LowSNR, *snr);
		return false;
	}

	*period = -1;

	amplitude->value = amax;

	if ( _streamConfig[_usedComponent].gain != 0.0 )
		amplitude->value /= _streamConfig[_usedComponent].gain;
	else {
		setStatus(MissingGain, 0.0);
		return false;
	}

	// - convert to micrometer
	amplitude->value *= 1E06;

	// - estimate peak-to-peak from absmax amplitude
	amplitude->value *= 2.0;

	return true;
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
void Detector::process(const Record *record, const DoubleArray &filteredData) {
	_amplProc.pickIndex = 0;

	SimpleDetector::process(record, filteredData);

	if ( _amplProc.isRunning() ) {
		calculateMaxAmplitude(record, _amplProc.pickIndex, filteredData.size(), filteredData);
		if ( _amplProc.isFinished() )
			sendMaxAmplitude(record);
	}
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor_Mjma::deconvolveData(Response *resp,
                                             DoubleArray &data,
                                             int numberOfIntegrations) {
	Math::Restitution::FFT::TransferFunctionPtr tf =
		resp->getTransferFunction(numberOfIntegrations);

	if ( tf == NULL ) {
		setStatus(DeconvolutionFailed, 0);
		return false;
	}

	
	Math::SeismometerResponse::Seismometer5sec paz(Math::Velocity);
	Math::Restitution::FFT::PolesAndZeros seis5sec(paz);

	Math::Restitution::FFT::TransferFunctionPtr cascade =
		*tf / seis5sec;

	// Remove linear trend
	double m,n;
	Math::Statistics::computeLinearTrend(data.size(), data.typedData(), m, n);
	Math::Statistics::detrend(data.size(), data.typedData(), m, n);

	return Math::Restitution::transformFFT(data.size(), data.typedData(),
	                                         _stream.fsamp, cascade.get(),
	                                         _config.respTaper, _config.respMinFreq, _config.respMaxFreq);
}

// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
bool AmplitudeProcessor_Md::deconvolveData(Response* resp,
                                           DoubleArray& data,
                                           int numberOfIntegrations) {
	if ( numberOfIntegrations < -1 )
		return false;

	SEISCOMP_DEBUG("Inside deconvolve function");

	double m, n;
	Math::Restitution::FFT::TransferFunctionPtr tf =
		resp->getTransferFunction(numberOfIntegrations < 0 ? 0 : numberOfIntegrations);

	if ( !tf )
		return false;

	Math::GroundMotion gm;

	if ( numberOfIntegrations < 0 )
		gm = Math::Displacement;
	else
		gm = Math::Velocity;

	Math::Restitution::FFT::TransferFunctionPtr cascade;
	Math::SeismometerResponse::WoodAnderson woodAndersonResp(gm, _config.woodAndersonResponse);
	Math::SeismometerResponse::Seismometer5sec seis5sResp(gm);
	Math::SeismometerResponse::L4C_1Hz l4c1hzResp(gm);

	Math::Restitution::FFT::PolesAndZeros woodAnderson(woodAndersonResp);
	Math::Restitution::FFT::PolesAndZeros seis5sec(seis5sResp);
	Math::Restitution::FFT::PolesAndZeros l4c1hz(l4c1hzResp);

	SEISCOMP_DEBUG("SEISMO = %d", aFile.SEISMO);

	switch ( aFile.SEISMO ) {
		case 1:
			cascade = *tf / woodAnderson;
		break;
		case 2:
			cascade = *tf / seis5sec;
		break;
		case 9:
			SEISCOMP_INFO("%s Applying filter L4C 1Hz to data", AMPTAG);
			cascade = *tf / l4c1hz;
		break;
		default:
			cascade = tf;
			SEISCOMP_INFO("%s No seismometer specified, no signal reconvolution performed", AMPTAG);
			return false;
		break;
	}

	// Remove linear trend
	Math::Statistics::computeLinearTrend(data.size(), data.typedData(), m, n);
	Math::Statistics::detrend(data.size(), data.typedData(), m, n);

	return Math::Restitution::transformFFT(data.size(), data.typedData(),
	    _stream.fsamp, cascade.get(), _config.respTaper, _config.respMinFreq,
	    _config.respMaxFreq);
}

int CreateCalibration( DoubleArray& Nexp, DoubleArray& teta, CalibrationParams& cal )
{
	CalibrationFuncParams in_params(Nexp.GetSize(), Nexp, teta, cal.n_p, cal.n_s, cal.alfa );
	CalibrationSolver FindFI( in_params );
	if( (cal.status=FindFI.Run(-45*DEGREE, 45*DEGREE, 1e-12))==GSL_SUCCESS ) 
	{
		double *A = FindFI.fparams.GetA(), *B = FindFI.fparams.GetB(), *N = FindFI.fparams.N;

		cal.fi0=FindFI.root;
		cal.L =	((N[1] - N[3])*(A[0] - A[3]) - (N[0] - N[3])*(A[1] - A[3])) / 
			((B[0] - B[3])*(A[1] - A[3]) - (B[1] - B[3])*(A[0] - A[3])); 
		cal.d0 = ((N[3] - N[0]) - cal.L*(B[0] - B[3]))/(A[0] - A[3]);
		cal.N0 = N[0] + cal.d0*A[0] + cal.L*B[0];

		cal.dt=FindFI.dt;
		cal.func_call_cntr=CalibrationSolver::func_call_cntr;
		cal.epsabs=FindFI.epsabs; cal.epsrel=FindFI.epsrel;
		cal.Nexp.RemoveAll(); cal.teta.RemoveAll();
		for(int i=0;i<Nexp.GetSize();i++) cal.Nexp.Add(N[i]);
		for(int i=0;i<teta.GetSize();i++) cal.teta.Add(teta[i]);
	}	
	return cal.status;
}