C++ MatrixXf::rows方法代码示例

    void operator()(
    const Eigen::MatrixXf &A,
    const Eigen::MatrixXf &B,
    Eigen::MatrixXf &Transform){
        assert(A.rows() == B.rows());
        assert(A.cols() == numColumnElements);
        assert(B.cols() == numColumnElements);

        Transform.resize(numColumnElements, numColumnElements);
        /*
        Transform.col(0) = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B.col(0));
        Transform.col(1) = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B.col(1));
        Transform.col(2) = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B.col(2));
        */
        Transform = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B);
    }

void
computeHistogram (const Eigen::MatrixXf &data, Eigen::MatrixXf &histogram, size_t bins, float min, float max)
{
    float bin_size = (max-min) / bins;
    int num_dim = data.cols();
    histogram = Eigen::MatrixXf::Zero (bins, num_dim);

    for (int dim = 0; dim < num_dim; dim++)
    {
        omp_lock_t bin_lock[bins];
        for(size_t pos=0; pos<bins; pos++)
            omp_init_lock(&bin_lock[pos]);

    #pragma omp parallel for firstprivate(min, max, bins) schedule(dynamic)
        for (int j = 0; j<data.rows(); j++)
        {
            int pos = std::floor( (data(j,dim) - min) / bin_size);

            if(pos < 0)
                pos = 0;

            if(pos > (int)bins)
                pos = bins - 1;

            omp_set_lock(&bin_lock[pos]);
            histogram(pos,dim)++;
            omp_unset_lock(&bin_lock[pos]);
        }

        for(size_t pos=0; pos<bins; pos++)
            omp_destroy_lock(&bin_lock[pos]);
    }
}

void coordinateDescentLasso(
        const Eigen::MatrixXf &data,
        const Eigen::VectorXf &output,
        Eigen::VectorXf &weights, 
        const float lambda,
        const int nIters,
        const bool verbose)
{
    const int nExamples = data.rows();
    const int nFeatures = data.cols();
    for(int iter = 0; iter < nIters; ++iter){
        const int featureInd = iter % nFeatures;
        float rho = 0;
        for(int i = 0; i < nExamples; ++i)
            rho += residualWithoutKWeight(
                    weights,
                    data.row(i).transpose(),
                    output[i],
                    featureInd) * data(i, featureInd);
        auto column = data.col(featureInd);
        float sumOfColumn = (column.transpose() * column).sum();
        weights[featureInd] = coordinateDescentStepLasso(weights[featureInd], sumOfColumn, rho, lambda);
        if(verbose){
            const float err = rss(weights, data, output);
            std::cout << "iter " << iter << " err " << err << std::endl;
            std::cout << weights << std::endl;
        }
    }
}

template <typename PointInT, typename PointOutT> void
pcl::IntensitySpinEstimation<PointInT, PointOutT>::computeIntensitySpinImage (
      const PointCloudIn &cloud, float radius, float sigma, 
      int k,
      const std::vector<int> &indices, 
      const std::vector<float> &squared_distances, 
      Eigen::MatrixXf &intensity_spin_image)
{
  // Determine the number of bins to use based on the size of intensity_spin_image
  int nr_distance_bins = static_cast<int> (intensity_spin_image.cols ());
  int nr_intensity_bins = static_cast<int> (intensity_spin_image.rows ());

  // Find the min and max intensity values in the given neighborhood
  float min_intensity = std::numeric_limits<float>::max ();
  float max_intensity = -std::numeric_limits<float>::max ();
  for (int idx = 0; idx < k; ++idx)
  {
    min_intensity = (std::min) (min_intensity, cloud.points[indices[idx]].intensity);
    max_intensity = (std::max) (max_intensity, cloud.points[indices[idx]].intensity);
  }

  float constant = 1.0f / (2.0f * sigma_ * sigma_);
  // Compute the intensity spin image
  intensity_spin_image.setZero ();
  for (int idx = 0; idx < k; ++idx)
  {
    // Normalize distance and intensity values to: 0.0 <= d,i < nr_distance_bins,nr_intensity_bins
    const float eps = std::numeric_limits<float>::epsilon ();
    float d = static_cast<float> (nr_distance_bins) * std::sqrt (squared_distances[idx]) / (radius + eps);
    float i = static_cast<float> (nr_intensity_bins) * 
              (cloud.points[indices[idx]].intensity - min_intensity) / (max_intensity - min_intensity + eps);

    if (sigma == 0)
    {
      // If sigma is zero, update the histogram with no smoothing kernel
      int d_idx = static_cast<int> (d);
      int i_idx = static_cast<int> (i);
      intensity_spin_image (i_idx, d_idx) += 1;
    }
    else
    {
      // Compute the bin indices that need to be updated (+/- 3 standard deviations)
      int d_idx_min = (std::max)(static_cast<int> (floor (d - 3*sigma)), 0);
      int d_idx_max = (std::min)(static_cast<int> (ceil  (d + 3*sigma)), nr_distance_bins - 1);
      int i_idx_min = (std::max)(static_cast<int> (floor (i - 3*sigma)), 0);
      int i_idx_max = (std::min)(static_cast<int> (ceil  (i + 3*sigma)), nr_intensity_bins - 1);
   
      // Update the appropriate bins of the histogram 
      for (int i_idx = i_idx_min; i_idx <= i_idx_max; ++i_idx)  
      {
        for (int d_idx = d_idx_min; d_idx <= d_idx_max; ++d_idx)
        {
          // Compute a "soft" update weight based on the distance between the point and the bin
          float w = expf (-powf (d - static_cast<float> (d_idx), 2.0f) * constant - powf (i - static_cast<float> (i_idx), 2.0f) * constant);
          intensity_spin_image (i_idx, d_idx) += w;
        }
      }
    }
  }
}

int main(int argc, char** argv)
{
	std::string p_in = "";
	std::string p_out = "out.tsv";

	namespace po = boost::program_options;
	po::options_description desc;
	desc.add_options()
		("help", "produce help message")
		("density", po::value(&p_in), "filename for input density")
		("out", po::value(&p_out), "filename for smoothed output density")
	;

	po::variables_map vm;
	po::store(po::parse_command_line(argc, argv, desc), vm);
	po::notify(vm);
	if(vm.count("help")) {
		std::cerr << desc << std::endl;
		return 1;
	}

	Eigen::MatrixXf rho = density::LoadDensity(p_in);
	std::cout << "Loaded density dim=" << rho.rows() << "x" << rho.cols() << ", sum=" << rho.sum() << "." << std::endl;

	Eigen::MatrixXf result;
	{
		boost::timer::auto_cpu_timer t;
		result = density::DensityAdaptiveSmooth(rho);
	}

	density::SaveDensity(p_out, result);
	std::cout << "Wrote result to file '" << p_out << "'." << std::endl;

	return 1;
}

	void addEigenMatrixRow( Eigen::MatrixXf &m )
	{
		Eigen::MatrixXf temp=m;
		m.resize(m.rows()+1,m.cols());
		m.setZero();
		m.block(0,0,temp.rows(),temp.cols())=temp;
	}

int run() {
    // store cooccurrence in Eigen sparse matrix object
    REDSVD::SMatrixXf A;
    const int ncontext = read_cooccurrence(c_cooc_file_name, A, verbose);

    // read U matrix from file
    Eigen::MatrixXf V;
    read_eigen_truncated_matrix(c_input_file_V_name, V, dim);
    // read S matrix from file
    Eigen::VectorXf S;
    read_eigen_vector(c_input_file_S_name, S, dim, 1.0-eig);

    // checking the dimensions
    if (V.rows() != ncontext){
        throw std::runtime_error("size mismatch between projection V matrix and the number of context words!!");
    }

    // starting projection
    if (verbose) fprintf(stderr, "Running the projection...");
    const double start = REDSVD::Util::getSec();
    Eigen::MatrixXf embeddings = A * V * S.asDiagonal().inverse();
    if (norm) embeddings.rowwise().normalize();
    if (verbose) fprintf(stderr, "done in %.2f.\n",REDSVD::Util::getSec() - start);

    // write out embeddings
    const char *c_output_name = get_full_path(c_cooc_dir_name, c_output_file_name);
    if (verbose) fprintf(stderr, "writing infered word embeddings in %s\n", c_output_name);
    write_eigen_matrix(c_output_name, embeddings);
    free((char*)c_output_name);

    return 0;
}

// draw all the lines between aPts and Pts that have a corr>threshold.
// if aPtInd is in the range of aPts, then draw only the lines that comes from aPts[aPtInd]
void drawCorrLines(PlotLines::Ptr lines, const vector<btVector3>& aPts, const vector<btVector3>& bPts, const Eigen::MatrixXf& corr, float threshold, int aPtInd) {
  vector<btVector3> linePoints;
  vector<btVector4> lineColors;
  float max_corr = corr.maxCoeff(); // color lines by corr, where corr has been mapped [threshold,max_corr] -> [0,1]
  for (int i=0; i<corr.rows(); ++i)
    for (int j=0; j<corr.cols(); ++j)
    	if (corr(i,j) > threshold) {
    		if (aPtInd<0 || aPtInd>=corr.rows() || aPtInd==i) {
					linePoints.push_back(aPts[i]);
					linePoints.push_back(bPts[j]);
					float color_factor = (corr(i,j)-threshold)/(max_corr-threshold); //basically, it ranges from 0 to 1
					lineColors.push_back(btVector4(color_factor, color_factor,0,1));
				}
    	}
  lines->setPoints(linePoints, lineColors);
}

// Send a matrix to MATLAB
IGL_INLINE void igl::mlsetmatrix(Engine** mlengine, std::string name, const Eigen::MatrixXf& M)
{
  if (*mlengine == 0)
    mlinit(mlengine);
  
  mxArray *A = mxCreateDoubleMatrix(M.rows(), M.cols(), mxREAL);
  double *pM = mxGetPr(A);
  
  int c = 0;
  for(int j=0; j<M.cols();++j)
    for(int i=0; i<M.rows();++i)
      pM[c++] = double(M(i,j));
  
  engPutVariable(*mlengine, name.c_str(), A);
  mxDestroyArray(A);
}

IGL_INLINE void igl::mlgetmatrix(Engine** mlengine, std::string name, Eigen::MatrixXf& M)
{
  if (*mlengine == 0)
    mlinit(mlengine);
  
  unsigned long m = 0;
  unsigned long n = 0;
  std::vector<double> t;
  
  mxArray *ary = engGetVariable(*mlengine, name.c_str());
  if (ary == NULL)
  {
    m = 0;
    n = 0;
    M = Eigen::MatrixXf(0,0);
  }
  else
  {
    m = mxGetM(ary);
    n = mxGetN(ary);
    M = Eigen::MatrixXf(m,n);
    
    double *pM = mxGetPr(ary);
    
    int c = 0;
    for(int j=0; j<M.cols();++j)
      for(int i=0; i<M.rows();++i)
        M(i,j) = pM[c++];
  }
  
  mxDestroyArray(ary);
}

Eigen::MatrixXf
readDescrFromFile(const std::string &file, int padding, int rowSize)
{

    // check if file exists
    boost::filesystem::path path = file;
    if ( ! (boost::filesystem::exists ( path ) && boost::filesystem::is_regular_file(path)) )
        throw std::runtime_error ("Given file path to read Matrix does not exist!");

    std::ifstream in (file.c_str (), std::ifstream::in);

    int bufferSize = 819200;
    //int bufferSize = rowSize * 10;

    char linebuf[bufferSize];



    Eigen::MatrixXf matrix;
    matrix.resize(0,rowSize);
    int j=0;
    while(in.getline (linebuf, bufferSize)){
        int start_s=clock();
        std::string line (linebuf);
        std::vector < std::string > strs_2;
        boost::split (strs_2, line, boost::is_any_of (" "));
        matrix.conservativeResize(matrix.rows()+1,rowSize);
        for (int i = 0; i < strs_2.size()-1; i++)
            matrix (j, i) = static_cast<float> (atof (strs_2[i].c_str ()));
        j++;
        int stop_s=clock();
        std::cout << "time: " << (stop_s-start_s)/double(CLOCKS_PER_SEC)*1000 << std::endl;
    }
    return matrix;
}

Eigen::MatrixXf SumMipMap(const Eigen::MatrixXf& img_big)
{
	// size of original image
	const unsigned int w_big = img_big.rows();
	const unsigned int h_big = img_big.cols();
	// size of reduced image
	const unsigned int w_sma = w_big / Q;
	const unsigned int h_sma = h_big / Q;
	// the computed mipmap will have 2^i size
	if(Q*w_sma != w_big || Q*h_sma != h_big) {
		throw std::runtime_error("ERROR: Q and size does not match in function SumMipMap!");
	}
	Eigen::MatrixXf img_small(w_sma, h_sma);
	for(unsigned int y=0; y<h_sma; ++y) {
		const unsigned int y_big = Q*y;
		for(unsigned int x=0; x<w_sma; ++x) {
			const unsigned int x_big = Q*x;
			float sum = 0.0f;
			for(unsigned int i=0; i<Q; ++i) {
				for(unsigned int j=0; j<Q; ++j) {
					sum += img_big(x_big+j, y_big+i);
				}
			}
			img_small(x, y) = sum;
		}
	}
	return img_small;
}

std::vector<Eigen::Vector2f> RectGrid(const Eigen::MatrixXf& density)
{
	const float width = static_cast<float>(density.rows());
	const float height = static_cast<float>(density.cols());
	const float numf = density.sum();
	const float d = std::sqrt(float(width*height) / numf);
	const unsigned int Nx = static_cast<unsigned int>(std::ceil(width / d));
	const unsigned int Ny = static_cast<unsigned int>(std::ceil(height / d));
	const float Dx = width / static_cast<float>(Nx);
	const float Dy = height / static_cast<float>(Ny);
	const float Hx = Dx/2.0f;
	const float Hy = Dy/2.0f;

	std::vector<Eigen::Vector2f> seeds;
	seeds.reserve(Nx*Ny);
	for(unsigned int iy=0; iy<Ny; iy++) {
		float y = Hy + Dy * static_cast<float>(iy);
		for(unsigned int ix=0; ix<Nx; ix++) {
			float x = Hx + Dx * static_cast<float>(ix);
			seeds.push_back(Eigen::Vector2f(x, y));
		}
	}

	return seeds;
}

void sumAndNormalize( Eigen::MatrixXf & out, const Eigen::MatrixXf & in, const Eigen::MatrixXf & Q ) {
	out.resize( in.rows(), in.cols() );
	for( int i=0; i<in.cols(); i++ ){
		Eigen::VectorXf b = in.col(i);
		Eigen::VectorXf q = Q.col(i);
		out.col(i) = b.array().sum()*q - b;
	}
}

void Eigen_to_G(const Eigen::MatrixXf& EigenM, Matrix* _GM)
{
  Matrix& GM = *_GM;
  GM.Dimension(EigenM.rows(), EigenM.cols());
  for(int i=0; i<GM.rows; i++)
    for(int j=0; j<GM.cols; j++)
      GM[i][j] = EigenM(i, j);
}