C++ math::Matrix类代码示例

void
AAKR::normalize(Math::Matrix& mean, Math::Matrix& std)
{
    // Resize mean and standard deviation variables.
    mean.resizeAndFill(1, sampleSize(), 0);
    std.resizeAndFill(1, sampleSize(), 0);

    // Compute mean.
    for (unsigned i = 0; i < sampleSize(); i++)
        mean(i) = sum(m_data.get(0, m_num_values - 1, i, i)) / m_num_values;

    // Compute standard deviation.
    for (unsigned j = 0; j < sampleSize(); j++)
    {
        double sum = 0;

        // Sum of the power of two difference
        // between the value and the mean.
        for (unsigned i = 0; i < m_num_values; i++)
            sum += std::pow(m_data(i, j) - mean(j), 2);

        // Standard deviation.
        std(j) = std::sqrt(sum / m_num_values);

        // Normalize each member of the data set.
        for (unsigned i = 0; i < m_num_values; i++)
        {
            if (std(j))
                m_norm(i, j) = (m_data(i, j) - mean(j)) / std(j);
            else
                m_norm(i, j) = 0;
        }
    }
}

TEST(MatrixTest, DetTest)
{
    const Math::Matrix mat1(
        {
            { -0.95880162984708284f,  0.24004047608997131f, -0.78172309932665407f, -0.11604124457222834f },
            { -0.36230592086261376f, -0.75778166876017261f,  0.33041059404631740f, -1.06001391941094836f },
            {  0.00260215210936187f,  1.27485610196385113f, -0.26149859846418033f, -0.59669701186364876f },
            {  0.36899429848485432f,  3.01720896813933104f,  2.10311476609438719f, -1.68627076626448269f }
        }
    );

    const float expectedDet1 = 4.07415413729671f;

    float ret1 = mat1.Det();
    EXPECT_TRUE(Math::IsEqual(ret1, expectedDet1, TEST_TOLERANCE));

    const Math::Matrix mat2(
        {
            { -1.0860073221346871f,  0.9150354098189495f, -0.2723201933559999f,  0.2922832160271507f },
            { -1.0248331304801788f, -2.5081237461125205f, -1.0277123574586633f, -0.2254690663329798f },
            { -1.4227635282899367f, -0.0403846809122684f,  0.9216148477171653f,  1.2517067488015878f },
            { -0.1160254467152022f,  0.8270675274393656f,  1.0327218739781614f, -0.3674886870220400f }
        }
    );

    const float expectedDet2 = -6.35122307880942f;

    float ret2 = mat2.Det();
    EXPECT_TRUE(Math::IsEqual(ret2, expectedDet2, TEST_TOLERANCE));
}

 DWI2QBI (const Math::Matrix<value_type>& FRT_SHT, Math::Matrix<value_type>& normalise_SHT, const DWI::Shells& shells) :
   FRT_SHT (FRT_SHT), 
   normalise_SHT (normalise_SHT), 
   shells (shells),
   dwi (FRT_SHT.columns()),
   qbi (FRT_SHT.rows()),
   amps (normalise ? normalise_SHT.rows() : 0) { }

void run ()
{

  DWI::Tractography::Properties properties;

  DWI::Tractography::Writer<> writer (argument.back(), properties);

  for (size_t n = 0; n < argument.size()-1; n++) {
    Math::Matrix<float> M;
    try {
      M.load (argument[n]);
      if (M.columns() != 3)
        throw Exception ("file \"" + argument[n] + "\" does not contain 3 columns - ignored");

      DWI::Tractography::Streamline<float> tck (M.rows());
      for (size_t i = 0; i < M.rows(); i++) {
        tck[i].set (M (i,0), M (i,1), M (i,2));
      }
      writer (tck);
      writer.total_count++;
    }
    catch (Exception) { }
  }

}

void CalibrationWnd::calcAdcI2Curr()
{
    // Make linear assumption, e.i.
    //I = a*adc + b*1;
    Math::Matrix<2> XtX;
    Math::Vector<2> XtY;
    int sz = adcI.size();
    for ( int i=0; i<2; i++ )
    {
        QVector<int> * a;
        if ( i==0 )
            a = &adcI;
        else
            a = 0;
        for ( int j=0; j<2; j++ )
        {
            QVector<int> * b;
            if ( j==0 )
                b = &adcI;
            else
                b = 0;
            qreal v = 0.0;
            for ( int k=0; k<sz; k++ )
            {
                qreal va = ( a ) ? a->at( k ) : 1.0;
                qreal vb = ( b ) ? b->at( k ) : 1.0;
                v += va * vb;
            }
            XtX[i][j] = v;
        }

        qreal v = 0.0;
        for ( int k=0; k<sz; k++ )
        {
            qreal va = ( a ) ? a->at( k ) : 1.0;
            qreal vy = curr.at( k );
            v += va * vy;
        }
        XtY[i] = v;
    }
    // A = (XtX)^-1 * XtY;
    Math::Matrix<2> invXtX;
    invXtX = XtX.inv();
    Math::Vector<2> A;
    for ( int i=0; i<2; i++ )
    {
        qreal v = 0.0;
        for ( int j=0; j<2; j++ )
        {
            v += invXtX[i][j] * XtY[j];
        }
        A[i] = v;
    }
    aAdcI = A[0];
    bAdcI = A[1];
}

  void writeAsciiMatrix(const std::string& fname, const Math::Matrix<T,P,S>& M,
                        const std::string& meta, const bool trans = false) {
    Math::Range start(0,0);
    Math::Range end(M.rows(), M.cols());

    std::ofstream ofs(fname.c_str());
    if (!ofs.is_open())
      throw(std::runtime_error("Cannot open " + fname + " for writing."));
    MatrixWriteImpl<T,P,S,internal::BasicMatrixFormatter<T> >::write(ofs, M, meta, start, end, trans);
  }

void
TestMatrix::runSubTest18(double& res, double& expected, std::string& subTestName)
{
    expected = 1;
    subTestName = "simple_symmetric_invert";

#ifdef COSMO_LAPACK
    Math::SymmetricMatrix<double> mat(2, 2);
    mat(0, 0) = 2;
    mat(1, 1) = 3;
    mat(1, 0) = 1;

    mat.writeIntoTextFile("test_files/matrix_test_18_original.txt");

    Math::SymmetricMatrix<double> invMat = mat;
    invMat.invert();

    invMat.writeIntoTextFile("test_files/matrix_test_18_inverse.txt");

    Math::Matrix<double> prod = mat;
    prod *= invMat;
    prod.writeIntoTextFile("test_files/matrix_test_18_product.txt");

    res = 1;
    for(int i = 0; i < prod.rows(); ++i)
    {
        for(int j = 0; j < prod.rows(); ++j)
        {
            if(i == j)
            {
                if(!Math::areEqual(prod(i, j), 1.0, 1e-5))
                {
                    output_screen("FAIL! Diagonal element " << i << " must be 1 but it is " << prod(i, j) << std::endl);
                    res = 0;
                }
            }
            else
            {
                if(!Math::areEqual(prod(i, j), 0.0, 1e-5))
                {
                    output_screen("FAIL! Non-diagonal element " << i << " " << j << " must be 0 but it is " << prod(i, j) << std::endl);
                    res = 0;
                }
            }
        }
    }
#else
    output_screen_clean("This test (below) is skipped because Cosmo++ has not been linked to lapack" << std::endl);
    res = 1;
#endif
}

void run()
{
  InputBufferType dwi_buffer (argument[0], Image::Stride::contiguous_along_axis (3));
  Math::Matrix<cost_value_type> grad = DWI::get_valid_DW_scheme<cost_value_type> (dwi_buffer);

  size_t dwi_axis = 3;
  while (dwi_buffer.dim (dwi_axis) < 2) ++dwi_axis;
  INFO ("assuming DW images are stored along axis " + str (dwi_axis));

  Math::Matrix<cost_value_type> bmatrix;
  DWI::grad2bmatrix (bmatrix, grad);

  Math::Matrix<cost_value_type> binv (bmatrix.columns(), bmatrix.rows());
  Math::pinv (binv, bmatrix);

  int method = 1;
  Options opt = get_options ("method");
  if (opt.size()) method = opt[0][0];

  opt = get_options ("regularisation");
  cost_value_type regularisation = 5000.0;
  if (opt.size()) regularisation = opt[0][0];

  opt = get_options ("mask");
  Ptr<MaskBufferType> mask_buffer;
  Ptr<MaskBufferType::voxel_type> mask_vox;
  if (opt.size()){
    mask_buffer = new MaskBufferType (opt[0][0]);
    Image::check_dimensions (*mask_buffer, dwi_buffer, 0, 3);
    mask_vox = new MaskBufferType::voxel_type (*mask_buffer);
  }


  Image::Header dt_header (dwi_buffer);
  dt_header.set_ndim (4);
  dt_header.dim (3) = 6;
  dt_header.datatype() = DataType::Float32;
  dt_header.DW_scheme() = grad;

  OutputBufferType dt_buffer (argument[1], dt_header);

  InputBufferType::voxel_type dwi_vox (dwi_buffer);
  OutputBufferType::voxel_type dt_vox (dt_buffer);

  Image::ThreadedLoop loop ("estimating tensor components...", dwi_vox, 1, 0, 3);
  Processor processor (dwi_vox, dt_vox, mask_vox, bmatrix, binv, method, regularisation, loop.inner_axes()[0], dwi_axis);

  loop.run_outer (processor);
}

void
AAKR::computeDistance(Math::Matrix query)
{
    if (query.rows() != 1)
        throw std::runtime_error("unable to compute distance: reference is not row vector.");

    if ((unsigned)query.columns() != sampleSize())
        throw std::runtime_error("unable to compute distance: sample size does not match.");

    m_distances.fill(0.0);

    // Fill distances vector.
    for (unsigned i = 0; i < m_num_values; i++)
    {
        Math::Matrix q = query - m_norm.row(i);
        m_distances(i) = std::sqrt(sum(q * transpose(q)));
    };
}

void Render(Gfx::CGLDevice *device, Gfx::CModelFile *modelFile)
{
    device->BeginScene();

    Math::Matrix persp;
    Math::LoadProjectionMatrix(persp, Math::PI / 4.0f, (800.0f) / (600.0f), 0.1f, 100.0f);
    device->SetTransform(Gfx::TRANSFORM_PROJECTION, persp);

    Math::Matrix id;
    id.LoadIdentity();
    device->SetTransform(Gfx::TRANSFORM_WORLD, id);

    Math::Matrix viewMat;
    Math::LoadTranslationMatrix(viewMat, TRANSLATION);
    Math::Matrix rot;
    Math::LoadRotationXZYMatrix(rot, ROTATION);
    viewMat = Math::MultiplyMatrices(viewMat, rot);
    device->SetTransform(Gfx::TRANSFORM_VIEW, viewMat);

    const std::vector<Gfx::ModelTriangle> &triangles = modelFile->GetTriangles();

    Gfx::VertexTex2 tri[3];

    for (int i = 0; i < static_cast<int>( triangles.size() ); ++i)
    {
        device->SetTexture(0, GetTexture(triangles[i].tex1Name));
        device->SetTexture(1, GetTexture(triangles[i].tex2Name));
        device->SetTextureEnabled(0, true);
        device->SetTextureEnabled(1, true);

        device->SetMaterial(triangles[i].material);

        tri[0] = triangles[i].p1;
        tri[1] = triangles[i].p2;
        tri[2] = triangles[i].p3;

        device->DrawPrimitive(Gfx::PRIMITIVE_TRIANGLES, tri, 3);
    }

    device->EndScene();
}

void TSmatrix::testInversion()
{
	Math::Matrix m = Math::Matrix::zeros(4,4);
	for (int i=0; i<3; i++)
		for (int j=0; j<3; j++)
			if (j<=i)
				m(i,j) = 1;
	m(3,3) = 1;
	Math::Matrix minv(4,4);
	Math::Matrix::invert(m, minv);
	

	Math::Matrix identity = Math::Matrix::eye(4);
	Math::Matrix minv2(4,4);
	m.linearSolve(minv2, identity);

	minv.sub(minv2);
	for (int i=0;i<4;i++)
		for (int j=0;j<4;j++)
			TEST_ASSERT(fabs(minv(i,j)) < 0.0001); 
}

void run () 
{
  Math::Matrix<value_type> directions = DWI::Directions::load_cartesian<value_type> (argument[0]);

  size_t num_permutations = 1e8;
  Options opt = get_options ("permutations");
  if (opt.size())
    num_permutations = opt[0][0];

  Shared eddy_shared (directions, num_permutations);
  Thread::run (Thread::multi (Processor (eddy_shared)), "eval thread");

  auto& signs = eddy_shared.get_best_signs();

  for (size_t n = 0; n < directions.rows(); ++n) 
    if (signs[n] < 0)
      directions.row(n) *= -1.0;

  bool cartesian = get_options("cartesian").size();
  DWI::Directions::save (directions, argument[1], cartesian);
}

void verify_matrix (Math::Matrix<float>& in, const node_t num_nodes)
{
  if (in.rows() != in.columns())
    throw Exception ("Connectome matrix is not square (" + str(in.rows()) + " x " + str(in.columns()) + ")");
  if (in.rows() != num_nodes)
    throw Exception ("Connectome matrix contains " + str(in.rows()) + " nodes; expected " + str(num_nodes));

  for (node_t row = 0; row != num_nodes; ++row) {
    for (node_t column = row+1; column != num_nodes; ++column) {

      const float lower_value = in (column, row);
      const float upper_value = in (row, column);

      if (upper_value && lower_value && (upper_value != lower_value))
        throw Exception ("Connectome matrix is not symmetrical");

      if (!upper_value && lower_value)
        in (row, column) = lower_value;

      in (column, row) = 0.0f;

  } }
}

    void save_bvecs_bvals (const Image::Header& header, const std::string& path)
    {

      std::string bvecs_path, bvals_path;
      if (path.size() >= 5 && path.substr (path.size() - 5, path.size()) == "bvecs") {
        bvecs_path = path;
        bvals_path = path.substr (0, path.size() - 5) + "bvals";
      } else if (path.size() >= 5 && path.substr (path.size() - 5, path.size()) == "bvals") {
        bvecs_path = path.substr (0, path.size() - 5) + "bvecs";
        bvals_path = path;
      } else {
        bvecs_path = path + "bvecs";
        bvals_path = path + "bvals";
      }

      const Math::Matrix<float>& grad (header.DW_scheme());
      Math::Matrix<float> G (grad.rows(), 3);

      // rotate vectors from scanner space to image space
      Math::Matrix<float> D (header.transform());
      Math::Permutation p (4);
      int signum;
      Math::LU::decomp (D, p, signum);
      Math::Matrix<float> image2scanner (4,4);
      Math::LU::inv (image2scanner, D, p);
      Math::Matrix<float> rotation = image2scanner.sub (0,3,0,3);
      Math::Matrix<float> grad_G = grad.sub (0, grad.rows(), 0, 3);
      Math::mult (G, float(0.0), float(1.0), CblasNoTrans, grad_G, CblasTrans, rotation);

      // deal with FSL requiring gradient directions to coincide with data strides
      // also transpose matrices in preparation for file output
      std::vector<size_t> order = Image::Stride::order (header, 0, 3);
      Math::Matrix<float> bvecs (3, grad.rows());
      Math::Matrix<float> bvals (1, grad.rows());
      for (size_t n = 0; n < G.rows(); ++n) {
        bvecs(0,n) = header.stride(order[0]) > 0 ? G(n,order[0]) : -G(n,order[0]);
        bvecs(1,n) = header.stride(order[1]) > 0 ? G(n,order[1]) : -G(n,order[1]);
        bvecs(2,n) = header.stride(order[2]) > 0 ? G(n,order[2]) : -G(n,order[2]);
        bvals(0,n) = grad(n,3);
      }

      bvecs.save (bvecs_path);
      bvals.save (bvals_path);

    }

void run () 
{
  try {
    Math::Matrix<value_type> directions = DWI::Directions::load_cartesian<value_type> (argument[0]);
    report (str(argument[0]), directions);
  }
  catch (Exception& E) {
    Math::Matrix<value_type> directions (str(argument[0]));
    DWI::normalise_grad (directions);
    if (directions.columns() < 3) 
      throw Exception ("unexpected matrix size for DW scheme \"" + str(argument[0]) + "\"");

    print (str(argument[0]) + " [ " + str(directions.rows()) + " volumes ]\n");
    DWI::Shells shells (directions);

    for (size_t n = 0; n < shells.count(); ++n) {
      Math::Matrix<value_type> subset (shells[n].count(), 3);
      for (size_t i = 0; i < subset.rows(); ++i)
        subset.row(i) = directions.row(shells[n].get_volumes()[i]).sub(0,3);
      report ("\nb = " + str(shells[n].get_mean()), subset);
    }
  }
}