test_vector.cxx

InsightToolkit-1.4.0(有大量的优化算法程序)

  }

  {
    vnl_vector<float> v(3); v[0]=1.f; v[1] = 2.f; v[2] = 3.f;
    vnl_matrix<float> m = outer_product(v, v);
    TEST("outer_product",
         (m(0,0)==1 && m(0,1)==2 && m(0,2)==3 &&
          m(1,0)==2 && m(1,1)==4 && m(1,2)==6 &&
          m(2,0)==3 && m(2,1)==6 && m(2,2)==9), true);
  }
  {
    vnl_vector<float> v(3); v[0]=1.f; v[1] = 2.f; v[2] = 3.f;
    TEST("vnl_vector<float> v (1.f, 2.f, 3.f)", v.size(), 3);
    v.x() = 1.f;
    v.y() = 2.f;
    v.z() = 3.f;
    TEST("v.set_x(1) and v[0]", v[0], 1);
    TEST("v.set_y(2) and v[1]", v[1], 2);
    TEST("v.set_z(3) and v[2]", v[2], 3);
    vnl_vector<float> v1(3, 0.f); v1[0]=1.f;
    vcl_cout << "v1 = " << v1 << vcl_endl;
    vnl_vector<float> v2(3, 0.f); v2[1]=1.f;
    vcl_cout << "v2 = " << v2 << vcl_endl;
    vnl_vector<float> v3(3, 0.f); v3[0]=-0.5f; v3[2]=0.5f;
    vcl_cout << "v3 = " << v3 << vcl_endl;
    vcl_cout << "v1 - v2 = " << v1 - v2 << vcl_endl;
    double ang = angle(v1,v2);
    vcl_cout << "angle(v1,v2) = " << ang << vcl_endl;
    ang *= 180*vnl_math::one_over_pi;
    vcl_cout << "angle(v1,v2) in degrees = " << ang << vcl_endl;
    vcl_cout << "v1.size()=" << v1.size() << "\n";
    vcl_cout << "v2.size()=" << v2.size() << "\n";
    vcl_cout << "vnl_cross_2d(v1,v2) = " << vnl_cross_2d(v1,v2) << vcl_endl;
    vcl_cout << "vnl_cross_3d(v1,v2) = " << vnl_cross_3d(v1,v2) << vcl_endl;
    TEST_NEAR("angle(v1,v2)", ang, 90.0, 1e-15);
    double ang2 = angle(v1,v3);
    vcl_cout << "angle(v1,v3) = " << ang << vcl_endl;
    ang2 *= 180*vnl_math::one_over_pi;
    vcl_cout << "angle(v1,v3) in degrees = " << ang2 << vcl_endl;
    TEST_NEAR("angle(v1,v3)", ang2, 135.0, 1e-6);
#if 0
    TEST("squared_distance_2d", squared_distance_2d(v1,v2), 2);
    TEST("squared_distance_3d", squared_distance_3d(v1,v2), 2);
#endif
  }

  {
    float vvalues [] = {1,0,0,0};
    vnl_vector<float> v (4,4,vvalues);
    v[0] = 1;
    v[1] = 0;
    v[2] = 0;
    v[3] = 0;
    TEST("v.squared_magnitude", v.squared_magnitude(), 1);
    TEST("v.magnitude", v.magnitude(), 1);
    TEST("v.normalize", (v1 = 3.0f * v, v1.normalize(), v1), v);
  }

  TEST("vnl_vector_ssd",
       vnl_vector_ssd(vnl_vector<float>(4, 0.0f), vnl_vector<float>(4, 1.0f)),
       4.0);
}

void vnl_vector_test_matrix()
{
  int mvalues[] = {1,2,3,
                   4,5,6};                      // product with matrices
  vnl_matrix<int> m(2,3,6, mvalues);

  int v2values [] = {1,0};
  int v3values [] = {1,0,0};
  vnl_vector<int> v, v2(2,2,v2values), v3(3,3,v3values);
  TEST("v.pre_multiply(m)",
       ((v = v3),
        (v.pre_multiply(m)),
        (v.size()==2 && v(0)==1 && v(1)==4)), true);
  TEST("v.post_multiply(m)",
       ((v = v2),
        (v.post_multiply(m)),
        (v.size()==3 && v(0)==1 && v(1)==2 && v(2)==3)), true);
  TEST("v*=m",
       ((v = v2),
        (v *= m),
        (v.size()==3 && v(0)==1 && v(1)==2 && v(2)==3)), true);
  TEST("v2*m",
       ((v = v2 * m),
        (v.size()==3 && v(0)==1 && v(1)==2 && v(2)==3)), true);
  TEST("m*v3",
       ((v = m * v3),
        (v.size()==2 && v(0)==1 && v(1)==4)), true);
}

void vnl_vector_test_conversion()
{
  int i, d;
  bool check;
  {
    // convert from a vnl_vector to a block array:
    int v1values [] = {1,2,3, 4,5,6, 7,8,9, 10,11,12};
    vnl_vector<int> v1 (12, 12, v1values);
    const int* data = v1.data_block();
    {
      check = true;
      for (d = 0; d < 12; d++)
        if (data[d] != d+1)
          check = false;
    }
    TEST("(const int*) m.data_block", check, true);

    typedef int block [12];
    const block& v2 = *((const block*) data);
    {
      check = true;
      for (i = 0; i < 12; i++)
        if (v1(i) != v2[i])
          check = false;
    }
    TEST("matrix(i)==block[i]", check, true);

    // convert from a block array to a vnl_vector:
    block b1;
    for (i = 0; i < 12; i++)
      b1[i] = i;
    data = ((const int*) b1);
    {
      check = true;
      for (d = 0; d < 12; d++)
        if (data[d] != d)
          check = false;
    }
    TEST("(const int*) block", check, true);
    vnl_vector<int> b2(data, 12);
    {
      check = true;
      for (i = 0; i < 12; i++)
        if (b1[i] != b2(i))
          check = false;
    }
    TEST("block[i]==matrix(i)", check, true);
  }
#if 0
  {
    // convert from a vnl_vector to a block array:
    vnl_vector<double> v1 (12, 12,
                          1.0,2.0,3.0, 4.0,5.0,6.0,
                          7.0,8.0,9.0, 10.0,11.0,12.0);
    const double* data = v1.data_block();
    {
      check = true;
      for (d = 0; d < 12; d++)
     if (data[d] != d+1)
       check = false;
    }
    TEST("(const double*) m.data_block", check, true);

    typedef double block [12];
    block& v2 = *((block*) data);
    {
      check = true;
      for (i = 0; i < 12; i++)
     if (v1(i) != v2[i])
       check = false;
    }
    TEST("matrix(i)==block[i]", check, true);

    // convert from a block array to a vnl_vector:
    block b1;
    for (i = 0; i < 12; i++)
      b1[i] = i;
    data = ((const double*) b1);             // & in ((const double*) &b1)
    {                                                // is not needed
      check = true;
      for (d = 0; d < 12; d++)
     if (data[d] != d)
       check = false;
    }
    TEST("(const double*) block", check, true);
    vnl_vector<double> b2(data, 12);
    {
      check = true;
    for (i = 0; i < 12; i++)
      if (b1[i] != b2(i))
     check = false;
    }
    TEST("block[i]==matrix(i)", check, true);
  }
#endif
}

#ifndef TIMING
#define TIMING 0
#endif

#if TIMING
#include <vul/vul_timer.h>
void vnl_vector_test_two_nrm2_timing(unsigned size, unsigned long num)
{
  vnl_vector<double> a(size);
  for (unsigned i= 0; i < size; i++)
    a(i) = i / size;

  double c=0;
  vul_timer t;
  for (unsigned i = 0; i <num;i++)
    c+= vnl_c_vector<double>::two_nrm2(a.begin(), size);
  double time = t.real();
  vcl_cout <<" Time for finding the two_nrm2 of " << size
           <<"-D vectors " << num << "times  = " << time / 1000.0 << "s.\n";
}

void vnl_vector_test_euclid_dist_sq_timing(unsigned size, unsigned long num)
{
  vnl_vector<double> a(size);
  vnl_vector<double> b(size);
  for (unsigned i= 0; i < size; i++)
  {
    a(i) = i / size;
    b(i) = i * i / size;
  }

  double c=0;
  vul_timer t;
  for (unsigned i = 0; i <num;i++)
    c+= vnl_c_vector<double>::euclid_dist_sq(a.begin(), b.begin(), size);
  double time = t.real();
  vcl_cout <<" Time for finding the euclid_dist_sq of " << size
           <<"-D vectors " << num << "times  = " << time / 1000.0 << "s.\n";
}

void vnl_vector_test_timing()
{
  vnl_vector_test_two_nrm2_timing(20000,20000ul);
  vnl_vector_test_two_nrm2_timing(1000,400000ul);
  vnl_vector_test_two_nrm2_timing(100,4000000ul);
  vnl_vector_test_two_nrm2_timing(4,100000000ul);
  vnl_vector_test_euclid_dist_sq_timing(40000,10000ul);
  vnl_vector_test_euclid_dist_sq_timing(1000,400000ul);
  vnl_vector_test_euclid_dist_sq_timing(100,4000000ul);
  vnl_vector_test_euclid_dist_sq_timing(4,100000000ul);
}
#endif

void vnl_vector_test_leak()           // use top4.1 to watch for memory.
{                                     // remember to kill process.
  while (true) {
    vnl_vector_test_int();
    vnl_vector_test_matrix();
    vnl_vector_test_conversion();
  }
}

#ifndef LEAK
#define LEAK 0
#endif

void test_vector()
{
  vnl_vector_test_int();
  vnl_vector_test_float();
  vnl_vector_test_matrix();
  vnl_vector_test_conversion();
#if TIMING
  vnl_vector_test_timing();
#endif
#if LEAK
  vnl_vector_test_leak();
#endif
}


TESTMAIN(test_vector);