C++ ASSERT_NEAR函数代码示例

TEST_F(testLHCNoiseFB, constructor2) {

    auto lhcnfb = new LHCNoiseFB(1.0, 0.1, 0.9, 100, false);

    auto params = std::string(TEST_FILES) +
                  "/LHCNoiseFB/constructor/test2/";
    f_vector_t v;

    util::read_vector_from_file(v, params + "g.txt");
    // util::dump(lhcnfb->fG, "fG\n");
    ASSERT_EQ(v.size(), lhcnfb->fG.size());

    auto epsilon = 1e-8;
    for (unsigned int i = 0; i < v.size(); ++i) {
        auto ref = v[i];
        auto real = lhcnfb->fG[i];

        ASSERT_NEAR(ref, real, epsilon * std::max(fabs(ref), fabs(real)))
            << "Testing of fG failed on i " << i << std::endl;
    }

    delete lhcnfb;
}

TEST_F(testLHCNoiseFB, track1)
{

    auto long_tracker = new RingAndRfSection();

    auto lhcnfb = new LHCNoiseFB(1.0, 0.1, 0.9, 1);

    f_vector_t res;
    for (uint i = 0; i < 100; ++i) {
        long_tracker->track();
        Context::Slice->track();
        lhcnfb->track();
        res.push_back(lhcnfb->fX);
    }

    auto params =
        std::string(TEST_FILES) + "/LHCNoiseFB/track/test1/";
    f_vector_t v;

    util::read_vector_from_file(v, params + "x.txt");
    // util::dump(lhcnfb->fG, "fG\n");

    ASSERT_EQ(v.size(), res.size());

    auto epsilon = 1e-8;
    for (unsigned int i = 0; i < v.size(); ++i) {
        auto ref = v[i];
        auto real = res[i];

        ASSERT_NEAR(ref, real, epsilon * std::max(fabs(ref), fabs(real)))
                << "Testing of fX failed on i " << i << std::endl;
    }

    delete lhcnfb;
    delete long_tracker;
}

TEST(WingComponentSegment5, GetSegmentIntersection_BUG)
{
    const char* filename = "TestData/CS_SegIntersectionBUG.xml";
    ReturnCode tixiRet;
    TiglReturnCode tiglRet;

    TiglCPACSConfigurationHandle tiglHandle = -1;
    TixiDocumentHandle tixiHandle = -1;

    tixiRet = tixiOpenDocument(filename, &tixiHandle);
    ASSERT_EQ (SUCCESS, tixiRet);
    tiglRet = tiglOpenCPACSConfiguration(tixiHandle, "", &tiglHandle);
    ASSERT_EQ(TIGL_SUCCESS, tiglRet);

    double xsi = 0;
    TiglBoolean hasWarning;
    tiglRet = tiglWingComponentSegmentGetSegmentIntersection(tiglHandle, "wing_Cseg", "wing_Seg2", 0.0589568, 1., 0.35, 1., 1., &xsi, &hasWarning);
    ASSERT_EQ(TIGL_SUCCESS, tiglRet);
    ASSERT_EQ(TIGL_TRUE, hasWarning);
    ASSERT_NEAR(1.0, xsi, 1e-2);

    ASSERT_EQ(TIGL_SUCCESS, tiglCloseCPACSConfiguration(tiglHandle));
    ASSERT_EQ(SUCCESS, tixiCloseDocument(tixiHandle));
}

	void TestDrawListRandom(int numEntries) {

		DrawList<MyData> list;
		std::vector<double> configured;
		std::vector<double> drawn;

		const unsigned int cnt = numEntries;
		const unsigned int numDraw = cnt * 4096;
		double probSum = 0;

		list.resize(cnt);
		drawn.resize(cnt);
		configured.resize(cnt);

		// fill
		for (unsigned int i = 0; i < cnt; ++i) {
			double rnd = double(rand()) / double(RAND_MAX);
			configured[i] = rnd;
			list.set(i, MyData(i,i), rnd);
			probSum += rnd;
		}

		// draw
		for (unsigned int i = 0; i < numDraw; ++i) {
			MyData& d = list.draw();
			drawn[d.x]++;
		}

		// compare
		for (unsigned int i = 0; i < cnt; ++i) {
			double a = (configured[i] / probSum);
			double b = (drawn[i] / numDraw);
			ASSERT_NEAR(a, b, a*0.50);			// allow 50% difference between cfg and drawn
		}

	}

TEST (LowPassFilterTest, DoesBothAtOnce) {
    // make two sine waves, one second long
    KeyFinder::AudioData a;
    a.setChannels(1);
    a.setFrameRate(frameRate);
    a.addToSampleCount(frameRate);
    for (unsigned int i = 0; i < frameRate; i++) {
        float sample = 0.0;
        sample += sine_wave(i, highFrequency, frameRate, magnitude); // high freq
        sample += sine_wave(i, lowFrequency, frameRate, magnitude); // low freq
        a.setSample(i, sample);
    }

    KeyFinder::LowPassFilter* lpf = new KeyFinder::LowPassFilter(filterOrder, frameRate, cornerFrequency, filterFFT);
    KeyFinder::Workspace w;
    lpf->filter(a, w);
    delete lpf;

    // test for lower wave only
    for (unsigned int i = 0; i < frameRate; i++) {
        float expected = sine_wave(i, lowFrequency, frameRate, magnitude);
        ASSERT_NEAR(expected, a.getSample(i), tolerance);
    }
}

TEST (AudioDataTest, DownsamplerResamplesSineWave) {
  unsigned int frameRate = 10000;
  unsigned int frames = frameRate * 4;
  float freq = 20;
  float magnitude = 32768.0;
  unsigned int factor = 5;

  KeyFinder::AudioData a;
  a.setChannels(1);
  a.setFrameRate(frameRate);
  a.addToSampleCount(frames);
  for (unsigned int i = 0; i < frames; i++)
    a.setSample(i, sine_wave(i, freq, frameRate, magnitude));

  a.downsample(factor);

  unsigned int newFrameRate = frameRate / factor;
  unsigned int newFrames = frames / factor;

  ASSERT_EQ(newFrameRate, a.getFrameRate());
  ASSERT_EQ(newFrames, a.getSampleCount());
  for (unsigned int i = 0; i < newFrames; i++)
    ASSERT_NEAR(sine_wave(i, freq, newFrameRate, magnitude), a.getSample(i), magnitude * 0.05);
}

TEST(CalcRPSTest, ReverseLeft) {
	double twist_msg[3] = {17, 0, -4.273};

	/* test linear_RPS calculations */

	double actual_left_linear_RPS = CalcLinearRPS( twist_msg[0] );
	double actual_right_linear_RPS = CalcLinearRPS( twist_msg[0] );

//	double expected_left_linear_RPS = twist_msg[0] / CIRCUMFERENCE;
	double expected_left_linear_RPS = (17/CIRCUMFERENCE);
	double expected_right_linear_RPS = (twist_msg[0]/CIRCUMFERENCE);

	ASSERT_NEAR(actual_left_linear_RPS, expected_left_linear_RPS,	0.01);
	ASSERT_NEAR(actual_right_linear_RPS, expected_right_linear_RPS, 0.01);

	/* test angular_rps calculations */

	double actual_left_angular_RPS = CalcAngularRPS( twist_msg[2], 'L' );
	double actual_right_angular_RPS = CalcAngularRPS( twist_msg[2], 'R' );

	double expected_left_angular_RPS = twist_msg[2]/(CalcWheelDirection(twist_msg[2], 'L')*DISTANCE_TO_CENTER*2*PI*WHEEL_RADIUS);

	double expected_right_angular_RPS = twist_msg[2] /
				(	CalcWheelDirection( twist_msg[2], 'R' )*
					DISTANCE_TO_CENTER*2*PI*WHEEL_RADIUS );

	ASSERT_NEAR(actual_left_angular_RPS, expected_left_angular_RPS,
							0.01);
	ASSERT_NEAR(actual_right_angular_RPS, expected_right_angular_RPS,
							0.01);

	/* test actual total rps */

	double expected_left_RPS = actual_left_linear_RPS + actual_left_angular_RPS;
	double expected_right_RPS = actual_right_linear_RPS + actual_right_angular_RPS;
	
	double actual_left_RPS = CalcRPS( twist_msg, 'L' );
	double actual_right_RPS = CalcRPS( twist_msg, 'R' );

	ASSERT_NEAR(actual_left_RPS, expected_left_RPS, 0.01);
	ASSERT_NEAR(actual_right_RPS, expected_right_RPS, 0.01);
}

TEST_F(WingComponentSegmentSimple, getPointInternal_accuracy)
{
    int compseg = 1;
    // now we have do use the internal interface as we currently have no public api for this
    tigl::CCPACSConfigurationManager & manager = tigl::CCPACSConfigurationManager::GetInstance();
    tigl::CCPACSConfiguration & config = manager.GetConfiguration(tiglHandle);
    tigl::CCPACSWing& wing = config.GetWing(1);
    tigl::CCPACSWingComponentSegment& segment = (tigl::CCPACSWingComponentSegment&) wing.GetComponentSegment(compseg);

    gp_Pnt point = segment.GetPoint(0.25, 0.5);
    ASSERT_NEAR(point.X(), 0.5, 1e-7);
    ASSERT_NEAR(point.Y(), 0.5, 1e-7);

    point = segment.GetPoint(0.5, 0.5);
    ASSERT_NEAR(point.X(), 0.5, 1e-7);
    ASSERT_NEAR(point.Y(), 1.0, 1e-7);

    point = segment.GetPoint(0.75, 0.5);
    ASSERT_NEAR(point.X(), 0.625, 1e-7);
    ASSERT_NEAR(point.Y(), 1.5, 1e-7);
}

    void checkInDense(std::vector<int>& amdRowPtr, std::vector<int>& amdColIndices, std::vector<T>& amdVals)
    {
        uBLAS::mapped_matrix<T> sparseDense(csrMatrixC.num_rows, csrMatrixC.num_cols, 0);
        uBLAS::mapped_matrix<T> boostDense(csrMatrixC.num_rows, csrMatrixC.num_cols, 0);

        // boost sparse_prod cancels out zeros and hence reports more accurately non-zeros
        // In clSPARSE, spGeMM produces more non-zeros, and considers some zeros as nonzeros.
        // Therefore converting to dense and verifying the output in  dense format
        // Convert CSR to Dense

        for (int i = 0; i < amdRowPtr.size() - 1; i++)
        {
            // i corresponds to row index
            for (int j = amdRowPtr[i]; j < amdRowPtr[i + 1]; j++)
                sparseDense(i, amdColIndices[j]) = amdVals[j];
        }

        for (int i = 0; i < this->C.index1_data().size() - 1; i++)
        {
            for (int j = this->C.index1_data()[i]; j < this->C.index1_data()[i + 1]; j++)
                boostDense(i, this->C.index2_data()[j]) = this->C.value_data()[j];
        }

        bool brelativeErrorFlag = false;
        bool babsErrorFlag = false;
        
        for (int i = 0; i < csrMatrixC.num_rows; i++)
        {
            for (int j = 0; j < csrMatrixC.num_cols; j++)
            {
                //ASSERT_EQ(boostDense(i, j), sparseDense(i, j));
#ifdef _DEBUG_SpMxSpM_
                ASSERT_NEAR(boostDense(i, j), sparseDense(i, j), SPGEMM_PREC_ERROR);
#else
                if (fabs(boostDense(i, j) - sparseDense(i, j)) > SPGEMM_PREC_ERROR)
                {
                    babsErrorFlag = true;
                    SCOPED_TRACE("Absolute Error Fail");
                    break;
                }
#endif
            }
        }
        // Relative Error
        for (int i = 0; i < csrMatrixC.num_rows; i++)
        {
            for (int j = 0; j < csrMatrixC.num_cols; j++)
            {
                float diff  = fabs(boostDense(i, j) - sparseDense(i, j));
                float ratio = diff / boostDense(i, j);
#ifdef _DEBUG_SpMxSpM_
                // ratio is less than or almost equal to SPGEMM_REL_ERROR
                EXPECT_PRED_FORMAT2(::testing::FloatLE, ratio, SPGEMM_REL_ERROR);
#else
                if (diff / boostDense(i, j) > SPGEMM_REL_ERROR)
                {
                    brelativeErrorFlag = true;
                    SCOPED_TRACE("Relative Error Fail");
                    break;
                }
#endif
            }
        }//
#ifndef _DEBUG_SpMxSpM_
        if (brelativeErrorFlag)
        {
            ASSERT_FALSE(babsErrorFlag);
        }

        if (babsErrorFlag)
        {
            ASSERT_FALSE(brelativeErrorFlag);
        }
#endif
    }// end

TEST(Solution, LinearElastic2D)
{
    try {
        MeshLib::IMesh *msh = MeshGenerator::generateStructuredRegularQuadMesh(2.0, 2, .0, .0, .0);
        GeoLib::Rectangle* _rec = new GeoLib::Rectangle(Point(0.0, 0.0, 0.0),  Point(2.0, 2.0, 0.0));
        Geo::PorousMedia pm;
        pm.hydraulic_conductivity = new NumLib::TXFunctionConstant(1.e-11);
        pm.porosity = new NumLib::TXFunctionConstant(0.2);
        Geo::Solid solid;
        solid.density = 3e+3;
        solid.poisson_ratio = 0.2;
        solid.Youngs_modulus = 10e+9;
        pm.solidphase = &solid;
        DiscreteSystem dis(msh);
        FemLib::LagrangeFeObjectContainer feObjects(msh);
        Geo::FemLinearElasticProblem* pDe = defineLinearElasticProblem(dis, *_rec, pm, &feObjects);
        TimeStepFunctionConstant tim(.0, 10.0, 10.0);
        pDe->setTimeSteppingFunction(tim);
        // options
        BaseLib::Options options;
        BaseLib::Options* op_lis = options.addSubGroup("LinearSolver");
        op_lis->addOption("solver_type", "CG");
        op_lis->addOption("precon_type", "NONE");
        op_lis->addOptionAsNum("error_tolerance", 1e-10);
        op_lis->addOptionAsNum("max_iteration_step", 500);
        BaseLib::Options* op_nl = options.addSubGroup("Nonlinear");
        op_nl->addOption("solver_type", "Linear");
        op_nl->addOptionAsNum("error_tolerance", 1e-6);
        op_nl->addOptionAsNum("max_iteration_step", 500);

        MyFunctionDisplacement f_u;
        f_u.define(&dis, pDe, &pm, options);
        f_u.setOutputParameterName(0, "u_x");
        f_u.setOutputParameterName(1, "u_y");
        f_u.setOutputParameterName(2, "Strain");
        f_u.setOutputParameterName(3, "Stress");
//        MyFunctionStressStrain f_sigma;
//        f_sigma.setOutputParameterName(0, "strain_xx");

        SerialStaggeredMethod method(1e-5, 100);
        AsyncPartitionedSystem apart1;
        apart1.setAlgorithm(method);
        apart1.resizeOutputParameter(4);
        apart1.setOutputParameterName(0, "u_x");
        apart1.setOutputParameterName(1, "u_y");
        apart1.setOutputParameterName(2, "Strain");
        apart1.setOutputParameterName(3, "Stress");
        apart1.addProblem(f_u);
        apart1.connectParameters();

        TimeSteppingController timestepping;
        timestepping.setTransientSystem(apart1);

        //const double epsilon = 1.e-3;
        timestepping.setBeginning(.0);
        timestepping.solve(tim.getEnd());

//        const MyNodalFunctionScalar* r_f_ux = apart1.getOutput<MyNodalFunctionScalar>(apart1.getOutputParameterID("u_x"));
//        const MyNodalFunctionScalar* r_f_uy = apart1.getOutput<MyNodalFunctionScalar>(apart1.getOutputParameterID("u_y"));
        const MyIntegrationPointFunctionVector* r_f_strain = apart1.getOutput<MyIntegrationPointFunctionVector>(apart1.getOutputParameterID("Strain"));
        const MyIntegrationPointFunctionVector* r_f_stress = apart1.getOutput<MyIntegrationPointFunctionVector>(apart1.getOutputParameterID("Stress"));
//        const IDiscreteVector<double>* vec_r_f_ux = r_f_ux->getDiscreteData();
//        const IDiscreteVector<double>* vec_r_f_uy = r_f_uy->getDiscreteData();
        const MyIntegrationPointFunctionVector::MyDiscreteVector* vec_strain = r_f_strain->getDiscreteData();
        const MyIntegrationPointFunctionVector::MyDiscreteVector* vec_stress = r_f_stress->getDiscreteData();

//        r_f_ux->printout();
//        r_f_uy->printout();
//        r_f_strain->printout();
//        r_f_stress->printout();

        const MathLib::LocalVector &strain1 = (*vec_strain)[0][0];
        const MathLib::LocalVector &stress1 = (*vec_stress)[0][0];
        double E = solid.Youngs_modulus;
        double nu = solid.poisson_ratio;
        double sx = .0; //E/((1.+nu)*(1-2*nu))*((1-nu)*ex+nu*ey);
        double sy = -1e+6; //E/((1.-nu)*(1-2*nu))*(nu*ex+(1-nu)*ey);
        double sxy = .0; //0.5*E/(1.+nu)*exy;
        double ex = (1+nu)/E*((1-nu)*sx-nu*sy);
        double ey = (1+nu)/E*((1-nu)*sy-nu*sx);
        double exy = 2*(1+nu)/E*sxy;
        double sz = nu*E/((1.+nu)*(1-2*nu))*(ex+ey);
        double epsilon = 1e-6;
        ASSERT_NEAR(ex, strain1(0), epsilon);
        ASSERT_NEAR(ey, strain1(1), epsilon);
        ASSERT_NEAR(.0, strain1(2), epsilon);
        ASSERT_NEAR(exy, strain1(3), epsilon);
        epsilon = 1;
        ASSERT_NEAR(sx, stress1(0), epsilon);
        ASSERT_NEAR(sy, stress1(1), epsilon);
        ASSERT_NEAR(sz, stress1(2), epsilon);
        ASSERT_NEAR(sxy, stress1(3), epsilon);

    } catch (const char* e) {
        std::cout << "***Exception caught! " << e << std::endl;
    }

}

TEST(macro_ASSERT_NEAR, floats_test_out_to_equal_within_given_epsilon)
{
  ASSERT_NEAR(3.1415, M_PI, 0.0001);
}

TEST_F(MeshLibProperties, AddDoublePointerProperties)
{
    ASSERT_TRUE(mesh != nullptr);
    std::string const& prop_name("GroupProperty");
    // check if a property with the name is already assigned to the mesh
    ASSERT_FALSE(mesh->getProperties().hasPropertyVector(prop_name));
    // data needed for the property
    const std::size_t n_prop_val_groups(10);
    const std::size_t n_items(mesh_size*mesh_size*mesh_size);
    std::vector<std::size_t> prop_item2group_mapping(n_items);
    // create simple mat_group to index mapping
    for (std::size_t j(0); j<n_prop_val_groups; j++) {
        std::size_t const lower(
            static_cast<std::size_t>(
                (static_cast<double>(j)/n_prop_val_groups)*n_items
            )
        );
        std::size_t const upper(
            static_cast<std::size_t>(
                (static_cast<double>(j+1)/n_prop_val_groups)*n_items
            )
        );
        for (std::size_t k(lower); k<upper; k++) {
            prop_item2group_mapping[k] = j;
        }
    }
    // obtain PropertyVector data structure
    boost::optional<MeshLib::PropertyVector<double*> &> group_properties(
        mesh->getProperties().createNewPropertyVector<double*>(
            prop_name, n_prop_val_groups, prop_item2group_mapping,
            MeshLib::MeshItemType::Cell
        )
    );
    ASSERT_EQ(prop_item2group_mapping.size(), (*group_properties).size());

    // initialize the property values
    for (std::size_t i(0); i<n_prop_val_groups; i++) {
        (*group_properties).initPropertyValue(i, static_cast<double>(i+1));
    }
    // check mapping to values
    for (std::size_t i(0); i<n_prop_val_groups; i++) {
        std::size_t const lower(
            static_cast<std::size_t>(
                (static_cast<double>(i)/n_prop_val_groups)*n_items
            )
        );
        std::size_t const upper(
            static_cast<std::size_t>(
                (static_cast<double>(i+1)/n_prop_val_groups)*n_items
            )
        );
        for (std::size_t k(lower); k<upper; k++) {
            ASSERT_NEAR(static_cast<double>(i+1), *(*group_properties)[k],
                std::numeric_limits<double>::epsilon());
        }
    }

    // the mesh should have the property assigned to cells
    ASSERT_TRUE(mesh->getProperties().hasPropertyVector(prop_name));
    // fetch the properties from the container
    boost::optional<MeshLib::PropertyVector<double*> const&>
        group_properties_cpy(mesh->getProperties().getPropertyVector<double*>(
            prop_name));
    ASSERT_FALSE(!group_properties_cpy);

    for (std::size_t k(0); k<n_items; k++) {
        ASSERT_EQ((*group_properties)[k], (*group_properties_cpy)[k]);
    }

    mesh->getProperties().removePropertyVector(prop_name);
    boost::optional<MeshLib::PropertyVector<double*> const&>
        removed_group_properties(mesh->getProperties().getPropertyVector<double*>(
            prop_name));

    ASSERT_TRUE(!removed_group_properties);
}

TEST(Math, VectorLength) {
	Vec3 v1(1,2,3);
	ASSERT_NEAR(3.7414, Math::length(v1), 0.001);
	Vec3 v2(0,0,-1);
	ASSERT_NEAR(1.0, Math::length(v2), 0.001);
}

TEST_F(NetworkControllerTest, trainNet){
	NetworkController nc;
	nc.loadNetwork(TEST_DATA_PATH("Student.na"));
	nc.loadNetwork(TEST_DATA_PATH("Student.sif"));
	nc.loadObservations(TEST_DATA_PATH("StudentData.txt"),TEST_DATA_PATH("controlStudent.json"));
	nc.trainNetwork();
	Network n = nc.getNetwork();
	Node grade = n.getNode("Grade");
	ASSERT_NEAR(0.3f,grade.getProbability(0,0),0.001);
	ASSERT_NEAR(0.05f,grade.getProbability(0,2),0.001);
	ASSERT_NEAR(0.9f,grade.getProbability(0,1),0.001);
	ASSERT_NEAR(0.5f,grade.getProbability(0,3),0.001);
	ASSERT_NEAR(0.4f,grade.getProbability(1,0),0.001);
	ASSERT_NEAR(0.25f,grade.getProbability(1,2),0.001);
	ASSERT_NEAR(0.08f,grade.getProbability(1,1),0.001);
	ASSERT_NEAR(0.3f,grade.getProbability(1,3),0.001);
	ASSERT_NEAR(0.3f,grade.getProbability(2,0),0.001);
	ASSERT_NEAR(0.7f,grade.getProbability(2,2),0.001);
	ASSERT_NEAR(0.02f,grade.getProbability(2,1),0.001);
	ASSERT_NEAR(0.2f,grade.getProbability(2,3),0.001);
	Node letter = n.getNode("Letter");
	ASSERT_NEAR(0.1f,letter.getProbability(0,0),0.001);
	ASSERT_NEAR(0.4f,letter.getProbability(0,1),0.001);
	ASSERT_NEAR(0.99f,letter.getProbability(0,2),0.001);
	ASSERT_NEAR(0.9f,letter.getProbability(1,0),0.001);
	ASSERT_NEAR(0.6f,letter.getProbability(1,1),0.001);
	ASSERT_NEAR(0.01f,letter.getProbability(1,2),0.001);
	Node intelligence = n.getNode("Intelligence");
	ASSERT_NEAR(0.7, intelligence.getProbability(0,0),0.001);
	ASSERT_NEAR(0.3, intelligence.getProbability(1,0),0.001);
	Node difficulty = n.getNode("Difficulty");
	ASSERT_NEAR(0.6, difficulty.getProbability(0,0),0.001);
	ASSERT_NEAR(0.4, difficulty.getProbability(1,0),0.001);	
	Node sat = n.getNode("SAT");
	ASSERT_NEAR(0.95f,sat.getProbability(0,0),0.001);
	ASSERT_NEAR(0.05f,sat.getProbability(1,0),0.001);
	ASSERT_NEAR(0.2f,sat.getProbability(0,1),0.001);
	ASSERT_NEAR(0.8f,sat.getProbability(1,1),0.001);
}

//.........这里部分代码省略.........

  rmd::SE3<float> T_world_ref;
  dataset.readCameraPose(T_world_ref, ref_entry);

  const auto curr_entry = dataset(curr_ind);
  cv::Mat curr_img;
  dataset.readImage(curr_img, curr_entry);
  cv::Mat curr_img_flt;
  curr_img.convertTo(curr_img_flt, CV_32F, 1.0f/255.0f);

  rmd::SE3<float> T_world_curr;
  dataset.readCameraPose(T_world_curr, curr_entry);

  const float min_scene_depth = 0.4f;
  const float max_scene_depth = 1.8f;

  rmd::SeedMatrix seeds(ref_img.cols, ref_img.rows, cam);

  StopWatchInterface * timer = NULL;
  sdkCreateTimer(&timer);
  sdkResetTimer(&timer);
  sdkStartTimer(&timer);

  seeds.setReferenceImage(
        reinterpret_cast<float*>(ref_img_flt.data),
        T_world_ref.inv(),
        min_scene_depth,
        max_scene_depth);

  sdkStopTimer(&timer);
  double t = sdkGetAverageTimerValue(&timer) / 1000.0;
  printf("setReference image CUDA execution time: %f seconds.\n", t);

  cv::Mat initial_depthmap(ref_img.rows, ref_img.cols, CV_32FC1);
  seeds.downloadDepthmap(reinterpret_cast<float*>(initial_depthmap.data));

  cv::Mat initial_sigma_sq(ref_img.rows, ref_img.cols, CV_32FC1);
  seeds.downloadSigmaSq(reinterpret_cast<float*>(initial_sigma_sq.data));

  cv::Mat initial_a(ref_img.rows, ref_img.cols, CV_32FC1);
  seeds.downloadA(reinterpret_cast<float*>(initial_a.data));

  cv::Mat initial_b(ref_img.rows, ref_img.cols, CV_32FC1);
  seeds.downloadB(reinterpret_cast<float*>(initial_b.data));

  const float avg_scene_depth = (min_scene_depth+max_scene_depth)/2.0f;
  const float max_scene_sigma_sq = (max_scene_depth - min_scene_depth) * (max_scene_depth - min_scene_depth) / 36.0f;
  for(size_t r=0; r<ref_img.rows; ++r)
  {
    for(size_t c=0; c<ref_img.cols; ++c)
    {
      ASSERT_FLOAT_EQ(avg_scene_depth, initial_depthmap.at<float>(r, c));
      ASSERT_FLOAT_EQ(max_scene_sigma_sq, initial_sigma_sq.at<float>(r, c));
      ASSERT_FLOAT_EQ(10.0f, initial_a.at<float>(r, c));
      ASSERT_FLOAT_EQ(10.0f, initial_b.at<float>(r, c));
    }
  }

  // Test initialization of NCC template statistics

  // CUDA computation
  cv::Mat cu_sum_templ(ref_img.rows, ref_img.cols, CV_32FC1);
  seeds.downloadSumTempl(reinterpret_cast<float*>(cu_sum_templ.data));
  cv::Mat cu_const_templ_denom(ref_img.rows, ref_img.cols, CV_32FC1);
  seeds.downloadConstTemplDenom(reinterpret_cast<float*>(cu_const_templ_denom.data));

  // Host computation
  cv::Mat ocv_sum_templ(ref_img.rows, ref_img.cols, CV_32FC1);
  cv::Mat ocv_const_templ_denom(ref_img.rows, ref_img.cols, CV_32FC1);

  const int side = seeds.getPatchSide();
  for(size_t y=side; y<ref_img.rows-side/2; ++y)
  {
    for(size_t x=side; x<ref_img.cols-side/2; ++x)
    {
      double sum_templ    = 0.0f;
      double sum_templ_sq = 0.0f;
      for(int patch_y=0; patch_y<side; ++patch_y)
      {
        for(int patch_x=0; patch_x<side; ++patch_x)
        {
          const double templ = (double) ref_img_flt.at<float>( y-side/2+patch_y, x-side/2+patch_x );
          sum_templ += templ;
          sum_templ_sq += templ*templ;
        }
      }
      ocv_sum_templ.at<float>(y, x) = (float) sum_templ;
      ocv_const_templ_denom.at<float>(y, x) = (float) ( ((double)(side*side))*sum_templ_sq - sum_templ*sum_templ );
    }
  }
  for(size_t r=side; r<ref_img.rows-side/2; ++r)
  {
    for(size_t c=side; c<ref_img.cols-side/2; ++c)
    {
      ASSERT_NEAR(ocv_sum_templ.at<float>(r, c), cu_sum_templ.at<float>(r, c), 0.00001f);
      ASSERT_NEAR(ocv_const_templ_denom.at<float>(r, c), cu_const_templ_denom.at<float>(r, c), 0.001f);
    }
  }

}