C++ Matrix::Dim方法代码示例

        /**
         * @brief    Adjoint transform (no constructor calls)
         *
         * @param  m    DWT to transform
         * @param  res  Reconstructed signal
         */
        inline
        void
        Adjoint      (const Matrix <T> & m, Matrix <T> & res)
        {

            assert (   m.Dim (0) == _sl1
                    && m.Dim (1) == _sl2
                    && (_dim == 2 || m.Dim (2) == _sl3)
                    && m.Dim () == res.Dim ());

            /* function pointer */
            (this ->* idpwt) (m, res);

        }

/**
 * @brief         Hann window
 * 
 * @param   size  Side lengths
 * @param   t     Scaling factor
 * @return        Window
 */
template <class T> inline Matrix< std::complex<T> >
hannwindow (const Matrix<size_t>& size, const T& t) {
	
	size_t dim = size.Dim(0);
	assert (dim > 1 && dim < 4);
	
	Matrix<double> res;
	
	if      (dim == 1) res = Matrix<double> (size[0], 1);
	else if (dim == 2) res = Matrix<double> (size[0], size[1]);
	else               res = Matrix<double> (size[0], size[1], size[2]);
	
	float          h, d;
	float          m[3];
	
	if (isvec(res)) {
		
		m[0] = 0.5 * size[0];
		m[1] = 0.0;
		m[2] = 0.0;
		
	} else if (is2d(res)) {
		
		m[0] = 0.5 * size[0];
		m[1] = 0.5 * size[1];
		m[2] = 0.0;
		
	} else {
		
		m[0] = 0.5 * size[0];
		m[1] = 0.5 * size[1];
		m[2] = 0.5 * size[2];
		
	}
	
	res = squeeze(res);
	
	for (size_t s = 0; s < res.Dim(2); s++)
		for (size_t r = 0; r < res.Dim(1); r++)
			for (size_t c = 0; c < res.Dim(0); c++) {
				d = pow( (float)pow(((float)c-m[0])/m[0],2) + pow(((float)r-m[1])/m[1],2) + pow(((float)s-m[2])/m[2],2) , (float)0.5);
				h = (d < 1) ? (0.5 + 0.5 * cos (PI * d)) : 0.0;
				res(c,r,s) = t * h;
			}
	
	return res;
	
}

	/**
	 * @brief    Backward transform
	 *
	 * @param  m To transform
	 * @return   Transform
	 */
	template <class T> inline static Matrix<T> 
	Adjoint    (const Matrix<T>& m) {

		size_t M = m.Dim(0);
		size_t N = m.Dim(1);
		
		Matrix<T> resx (M, N);
		Matrix<T> resy (M, N);

#pragma omp for
			for (size_t i = 0; i < M; i++) {
				resy (0,i) = -m(0,i,0);
				for (size_t j = 1; j < N-1; j++)
					resy (j,i) = m(j-1,i,0) - m(j,i,0);
				resy (N-1,i) = m(N-2,i,0);
			}
			
#pragma omp for
			for (size_t i = 0; i < N; i++) {
				resx (i,0) = -m(i,0,1);
				for (size_t j = 1; j < M-1; j++)
					resx (i,j) = m(i,j-1,1) - m(i,j,1);
				resx (i,M-1) = m(i,M-2,1);
			}
		
		resx += resy;
		
		return resx;
		
	}

	/**
	 * @brief    2D Forward transform
	 *
	 * @param  m To transform
	 * @return   Transform
	 */
	template <class T> inline static Matrix<T> 
	Trafo (const Matrix<T>& m) {
		
		size_t M = m.Dim(0);
		size_t N = m.Dim(1);
		
		Matrix<T> res (M, N, 2); 
		
#pragma omp for
			for (size_t i = 0; i < N; i++) {
				for (size_t j = 0; j < M-1; j++)
					res(j,i,0) = m(j+1,i) - m(j,i);
				res (M-1,i,0) = T(0.0);
			}
			
#pragma omp for
			for (size_t i = 0; i < M; i++) {
				for (size_t j = 0; j < N-1; j++)
					res(i,j,1) = m(i,j+1) - m(i,j);
				res (i,N-1,1) = T(0.0);
			}
		
		return res;
		
	}	

 template<class T> inline static Matrix<T>
 exp (const Matrix<T>& m) {
     Matrix<T> ret (m.Dim(), m.Res());
     #pragma omp parallel for
     for (int i = 0; i < numel(m); i++)
         ret[i] = std::exp(m[i]);
     return ret;
 }

/**
 * @brief       Which elements are Inf
 *
 * @param  M    Matrix
 * @return      Matrix of booleans true where inf
 */
template <class T> inline static  Matrix<cbool>
isnan (const Matrix<T>& M) {

    Matrix<cbool> res (M.Dim());
    for (size_t i = 0; i < res.Size(); ++i)
        res.Container()[i] = (is_nan(TypeTraits<T>::Imag(M[i]))||is_nan(TypeTraits<T>::Imag(M[i])));
    return res;

}

/**
 * @brief     Is matrix X-dimensional?
 *
 * @param  M  Matrix
 * @param  d  Dimension
 * @return    X-dimensional?
 */
template <class T>  inline static  bool isxd (const Matrix<T>& M, size_t d) {

    size_t l = 0;

    for (size_t i = 0; i < M.NDim(); ++i)
        if (M.Dim(i) > 1) ++l;

    return (l == d);

}

/**
 * @brief      Mean reducing a dimension
 *
 * @param  M   Matrix
 * @param  d   Dimension
 * @return     Average of M reducing d matrix
 */
template <class T> static inline Matrix<T>
mean (const Matrix<T>& M, size_t d) {
	
	Matrix<T> res  = M;
	float     quot = (float) res.Dim(d);
		
	res = sum (res, d);
	
	return res / quot;
	
}

/**
 * @brief          FLip left right
 *
 * @param  M        Matrix
 * @return         Flipped matrix
 */
template <class T> inline static Matrix<T> fliplr (const Matrix<T>& M)  {

    size_t srow = size(M,1), scol = size(M,0), nrow = numel (M)/srow;
    Matrix<T> res (M.Dim());

    for (size_t i = 0; i < nrow; ++i)
        for (size_t j = 0; j < srow; ++j)
            res[j*scol+i] = M[(srow-1-j)*scol+i];

    return res;

}

void Matrix::MatrixConcat(Matrix& concat_mat)
{
	if (this->NumElements() == 0)
	{
		assert(0);
	/*	this->Initialize(concat_mat.Dims());
		for (int i = 0; i < concat_mat.NumElements(); ++i)
		{
			this->tensor_array->at(i) = concat_mat.tensor_array->at(i);
		}*/

	}
	else
	{
		assert(concat_mat.Dim(1) == this->Dim(1)); // assert they have same number of columns
		for (int i = 0; i  < concat_mat.NumElements(); ++i)
		{
			this->tensor_array->push_back(concat_mat.tensor_array->at(i));
		}
		this->dims->at(0) += concat_mat.Dim(0);
	}
}

/**
 * @brief       nxnxn cube with sphere centered at p
 *
 * @param  p    Center point of sphere
 * @param  n    Side length of cube
 * @param  s    Scaling factor
 * @return      Matrix with circle
 */
template <class T> inline static Matrix<T>
sphere (const float* p, const size_t n, const T s = T(1)) {

	Matrix<T> res (n,n,n);

	float m[3];
	float rad;

	rad = p[0] * float(n) / 2.0;

	m[0] = (1.0 - p[1]) * float(n) / 2.0;
	m[1] = (1.0 - p[2]) * float(n) / 2.0;
	m[2] = (1.0 - p[3]) * float(n) / 2.0;

	for (size_t s = 0; s < res.Dim(2); s++)
		for (size_t r = 0; r < res.Dim(1); r++)
			for (size_t c = 0; c < res.Dim(0); c++)
				res(c,r) = ( pow (((float)c-m[0])/rad, (float)2.0) + pow (((float)r-m[1])/rad, (float)2.0) + pow (((float)s-m[2])/rad, (float)2.0) <= 1.0) ? s : T(0.0);

	return res;

}

		template <class T> bool
		Write (const Matrix<T>& M, const std::string& uri = "") {

			assert (m_file != NULL);

			dtype dt = CODTraits<T>::dt;
			size_t n, l;

			// Dump type
			if (!mwrite(&dt,         1, m_file, "data type"))
				return false;

			n = M.NDim();
			if (!mwrite(&n,          1, m_file, "data dimensions"))
				return false;

			// Dump dimensions
			if (!mwrite(M.Dim(),        m_file, "dimensions"))
				return false;

			// Dump resolutions
			if (!mwrite(M.Res(),        m_file, "resolutions"))
				return false;

			// Size of name and name
			n = uri.size();
			if (!mwrite(&n,          1, m_file, "name length"))
				return false;

			// Dump name
			if (!mwrite(uri.c_str(), n, m_file, "name"))
				return false;

			// Dump data
			if (!mwrite(M.Container(),  m_file, "data"))
				return false;

			if (!mwrite(delim.c_str(), delim.length(), m_file, "delimiter"))
				return false;

			return true;

		}

/**
 * @brief      3D resample data to new grid size
 *
 * @param  M   Incoming data
 * @param  f   Resampling factor in all 3 dimensions
 * @param  im  Interpolation method (LINEAR|BSPLINE)
 *
 * @return     Resampled data
 */
template<class T> static Matrix<T> 
resample (const Matrix<T>& M, const Matrix<double>& f, const InterpMethod& im) {
	

	Matrix <T> res = M;
	
#ifdef HAVE_INSIGHT
	
	typedef typename itk::OrientedImage< T, 3 > InputImageType;
	typedef typename itk::OrientedImage< T, 3 > OutputImageType;
	typedef typename itk::IdentityTransform< double, 3 > TransformType;
	typedef typename itk::LinearInterpolateImageFunction< InputImageType, double > InterpolatorType;
	typedef typename itk::ResampleImageFilter< InputImageType, InputImageType > ResampleFilterType;
	
	typename InterpolatorType::Pointer linterp = InterpolatorType::New();
	
	TransformType::Pointer trafo = TransformType::New();
	trafo->SetIdentity();
	
	typename InputImageType::SpacingType space;
	space[0] = 1.0/f[0];
	space[1] = 1.0/f[1];
	space[2] = 1.0/f[2];
	
	typedef typename InputImageType::SizeType::SizeValueType SizeValueType;
	typename InputImageType::SizeType size; 
	size[0] = static_cast<SizeValueType>(res.Dim(0));
	size[1] = static_cast<SizeValueType>(res.Dim(1));
	size[2] = static_cast<SizeValueType>(res.Dim(2));
	
	typename itk::OrientedImage< T, 3 >::Pointer input = itk::OrientedImage< T, 3 >::New();
	typename itk::OrientedImage< T, 3 >::Pointer output = itk::OrientedImage< T, 3 >::New();
	
	typename itk::Image< T, 3 >::IndexType ipos;
	ipos[0] = 0; ipos[1] = 0; ipos[2] = 0;
	typename itk::Image< T, 3 >::IndexType opos;
	opos[0] = 0; opos[1] = 0; opos[2] = 0;
	
	typename itk::Image< T, 3 >::RegionType ireg;
	ireg.SetSize(size);
	ireg.SetIndex(ipos);
	input->SetRegions(ireg);
	input->Allocate();
	
	typename itk::Image< T, 3 >::RegionType oreg;
	oreg.SetSize(size);
	ireg.SetIndex(opos);
	output->SetRegions(oreg);
	output->Allocate();
	
	for (size_t z = 0; z < res.Dim(2); z++)
		for (size_t y = 0; y < res.Dim(1); y++)
			for (size_t x = 0; x < res.Dim(0); x++) {
				ipos[0] = x; ipos[1] = y; ipos[2] = z;
				input->SetPixel (ipos, res.At(x,y,z));
			}
	
	typename ResampleFilterType::Pointer rs = ResampleFilterType::New();
	rs->SetInput( input );
	rs->SetTransform( trafo );
	rs->SetInterpolator( linterp );
	rs->SetOutputOrigin ( input->GetOrigin());
	rs->SetOutputSpacing ( space );
	rs->SetOutputDirection ( input->GetDirection());
	rs->SetSize ( size );
	rs->Update ();
	
	output = rs->GetOutput();
	
	res = Matrix<T> (res.Dim(0)*f[0], res.Dim(1)*f[1], res.Dim(2)*f[2]);
	res.Res(0) = res.Res(0)/f[0];
	res.Res(1) = res.Dim(1)/f[1];
	res.Res(2) = res.Dim(2)/f[2];
	
	for (size_t z = 0; z < res.Dim(2); z++)
		for (size_t y = 0; y < res.Dim(1); y++)
			for (size_t x = 0; x < res.Dim(0); x++) {
				opos[0] = x; opos[1] = y; opos[2] = z;
				res.At(x,y,z) = output->GetPixel (opos);
			}
	
#else 
	
	printf ("ITK ERROR - Resampling not performed without ITK!\n");
	
#endif
	
	return res;
	
}

		template<class T> bool
		Write (const Matrix<T>& M, const std::string& uri) {

			T t;
			Group group, *tmp;
			std::string path;

            boost::tokenizer<> tok(uri);
			std::vector<std::string> sv (Split (uri, "/"));
			std::string name = sv[sv.size() - 1];
			sv.pop_back(); // data name not part of path

			if (sv.size() == 0)
				path = "/";
			else
				for (size_t i = 0; i < sv.size(); i++) {
					if (sv[i].compare(""))
						path += "/";
						path += sv[i];
				}

			if (this->m_verb)
				printf ("Creating dataset %s at path (%s)\n", name.c_str(), path.c_str());

			try {

				group = m_file.openGroup(path);
				if (this->m_verb)
					printf ("Group %s opened for writing\n", path.c_str()) ;

			} catch (const Exception& e) {

				for (size_t i = 0, depth = 0; i < sv.size(); i++) {

					if (sv[i].compare("")) {

						try {
							group = (depth) ? (*tmp).openGroup(sv[i])   : m_file.openGroup(sv[i]);
						} catch (const Exception& e) {
							group = (depth) ? (*tmp).createGroup(sv[i]) : m_file.createGroup(sv[i]);
						}

						tmp = &group;
						depth++;

					}

				}

			}

			// One more field for complex numbers
			size_t tmpdim = ndims(M);

			std::vector<hsize_t> dims (tmpdim);

			for (size_t i = 0; i < tmpdim; i++)
				dims[i] = M.Dim(tmpdim-1-i);

			if (is_complex(t)) {
				dims.push_back(2);
				tmpdim++;
			}

			DataSpace space (tmpdim, &dims[0]);
			PredType*  type = HDF5Traits<T>::PType();

			DataSet set = group.createDataSet(name, (*type), space);

			set.write   (M.Ptr(), (*type));
			set.close   ();
			space.close ();

            return true;

		}

// default constructor
TensorJTree::TensorJTree(BayesianNetwork* bNet, Matrix& tree_matrix, vector<vector<int>*>& nodes_to_vars, int max_obs_per_mode, vector<int>& template_vec, int max_runs)
{
	elimination_order = NULL;
	this->bNet = bNet;
	this->tree_matrix = new Matrix(tree_matrix);
	this->all_obs_vars = new vector<int>();

	for (int i = 0; i < bNet->NumNodes(); ++i)
	{
		if (bNet->is_observed(i))
		{
			this->all_obs_vars->push_back(i);
		}
	}


	this->node_list = new vector<TensorJTNode*>();
	for (int i = 0; i < tree_matrix.Dim(0); ++i)
	{
		node_list->push_back(new TensorJTNode(this, i, max_obs_per_mode, max_runs));
	}

	// store parents and children

	for (int i = 0; i < NumCliques(); ++i)
	{
		vector<int> parents;
		vector<int> children;

		this->tree_matrix->RowFind(parents, i);
		this->tree_matrix->ColFind(children, i);

		if (parents.size() == 0)
			this->root = i;
		else
			node_list->at(i)->SetParent(parents[0]);

		node_list->at(i)->SetChildren(children);


	}

	// store C, S, and R vars
	for (int i = 0; i < NumCliques(); ++i)
	{
		vector<int> i_Sep;
		vector<int>& i_vars = *(nodes_to_vars[i]);

		int parent = GetParent(i);

		if (parent >= 0)
		{
			vector<int>& parent_vars = *(nodes_to_vars[parent]);
			VectorPlus::Intersect(i_Sep, i_vars, parent_vars);
		}

		node_list->at(i)->SetSVars(i_Sep);
	}

	for (int i = 0; i < NumCliques(); ++i)
	{
		vector<int> reduced_C_vars;
		VectorPlus::Union(reduced_C_vars, GetSVars(i));
		
		for (int j = 0; j < nodes_to_vars[i]->size(); ++j)
		{
			if (bNet->is_observed(nodes_to_vars[i]->at(j)))
			{
				if (!VectorPlus::Contains(reduced_C_vars, nodes_to_vars[i]->at(j)))
				{
					reduced_C_vars.push_back(nodes_to_vars[i]->at(j));
				}
			}
		}

		vector<int>& children = GetChildren(i);
		for (int j = 0; j < children.size(); ++j)
		{
			VectorPlus::Union(reduced_C_vars, GetSVars(children[j]));
		}
		node_list->at(i)->SetCVars(reduced_C_vars);
		
		int parent = GetParent(i);
		if (parent >= 0)
		{
			vector<int>& parent_vars = *(nodes_to_vars[parent]);
			vector<int> i_Res;
			VectorPlus::SetDiff(i_Res, reduced_C_vars, GetSVars(i));
			node_list->at(i)->SetRVars(i_Res);
		}
		else
		{
			node_list->at(i)->SetRVars(reduced_C_vars);
		}

	}
	
	for (int i = 0; i < NumCliques(); ++i)
	{
//.........这里部分代码省略.........