C++ Mat::col方法代码示例

bool gram_schmidt(Mat &src)
{
	Mat a(3, 1, CV_64FC1);
	Mat b(3, 1, CV_64FC1);
	Mat c(3, 1, CV_64FC1);


	src.col(0).copyTo(a);
	src.col(1).copyTo(b);
	src.col(2).copyTo(c);



	b = b - ((a.dot(b)) / (a.dot(a)))*a;



	c = c - a.dot(c) / a.dot(a)*a - b.dot(c)/b.dot(b)*b;




	norm_vec(a);
	norm_vec(b);
	norm_vec(c);

	a.copyTo(src.col(0));
	b.copyTo(src.col(1));
	c.copyTo(src.col(2));

	return true;

}

	void findBlobs::sortTheDistFromThreePoint(const Mat& ip, const Mat& mc, vector<int>& idx)
	{
		idx.clear();
		int numOfCol = mc.cols;

		Mat ip1 = ip.col(0);
		Mat ip1Mat = ip1 * Mat::ones(1,numOfCol, ip.type());
		subtract(ip1Mat, mc, ip1Mat);
		pow(ip1Mat, 2, ip1Mat);

		Mat ip2 = ip.col(1);
		Mat ip2Mat = ip2 * Mat::ones(1,numOfCol, ip.type());
		subtract(ip2Mat, mc, ip2Mat);
		pow(ip2Mat, 2, ip2Mat);

		Mat ip3 = ip.col(2);
		Mat ip3Mat = ip3 * Mat::ones(1,numOfCol, ip.type());
		subtract(ip3Mat, mc, ip3Mat);
		pow(ip3Mat, 2, ip3Mat);

		std::vector<double> dist;
		dist.clear();
		for (int i = 0; i < numOfCol; i++)
		{
			Scalar s = sum(ip1Mat.col(i)) + sum(ip2Mat.col(i)) + sum(ip3Mat.col(i));
			dist.push_back(s[0]);
			idx.push_back(i);
		}
		sort(idx.begin(), idx.end(),
			[& dist](size_t i1, size_t i2) {return dist[i1] <  dist[i2];});
	}

void sharpen(const Mat&image, Mat&result)
{
	//分配图像
	result.create(image.rows, image.cols, image.type());
	//处理除了第一行和最后一行
	for (int i = 1; i < image.rows - 1; i++){
		const uchar* previous =
			image.ptr<const uchar>(i - 1);
		const uchar* current =
			image.ptr<const uchar>(i);
		const uchar* next =
			image.ptr<const uchar>(i - 1);
		//output
		uchar*output = result.ptr<uchar>(i);
		for (int j = 1; j < image.cols - 1; j++){
			//saturate对结果进行截断
			*output++ = cv::saturate_cast<uchar>(
				5 * current[j] - current[j-1]
				- current[j + 1] - previous[j] - next[j]
				);
		}
	}
	//未处理的像素设置为0
	result.row(0).setTo(Scalar(0));
	result.row(result.rows - 1).setTo(Scalar(0));
	result.col(0).setTo(Scalar(0));

	result.col(result.cols - 1).setTo(Scalar(0));

}

Mat computeConvolution(Mat &m, double sigma, bool highpass) {

    // TODO only allow C1 or C3  CV_Assert(src.channels() == 3);
    Mat result;

    if (!m.empty()) {
        result = m.clone();
        // Store type to restore it when the convolution is computed
        int type = result.channels() == 1 ? CV_64F : CV_64FC3;
        // Convert the image to a 64F type, (one or three channels)
        result.convertTo(result, type);
        Mat kernel = myGetGaussianKernel1D(sigma, highpass);
        // This kernel is separable, apply convolution for rows and columns separately
        for (int i = 0; i < result.rows; i++) {
            Mat row = result.row(i);
            row = convolutionOperator1D(row, kernel, BORDER_REFLECT);
            row.copyTo(result.row(i));
        }
        for (int i = 0; i < result.cols; i++) {
            Mat col = result.col(i);
            col = convolutionOperator1D(col, kernel, BORDER_REFLECT);
            col.copyTo(result.col(i));
        }
        result.convertTo(result, m.type());
    }
    return result;
}

/*step的方式遍历*/
void sharpen2(const Mat &image, Mat &result) 
{
	result.create(image.size(), image.type()); // allocate if necessary

	int step= image.step1();
	const uchar* previous= image.data;		// ptr to previous row
	const uchar* current=  image.data+step; // ptr to current row
	const uchar* next= image.data+2*step;   // ptr to next row
	uchar *output= result.data+step;		// ptr to output row

	for (int j= 1; j<image.rows-1; j++) { // for each row (except first and last)
		for (int i=1; i<image.cols-1; i++) { // for each column (except first and last)

			output[i]= saturate_cast<uchar>(5*current[i]-current[i-1]-current[i+1]-previous[i]-next[i]); 
		}

		previous+= step;
		current+= step;
		next+= step;
		output+= step;
	}

	// Set the unprocess pixels to 0
	result.row(0).setTo(cv::Scalar(0));
	result.row(result.rows-1).setTo(cv::Scalar(0));
	result.col(0).setTo(cv::Scalar(0));
	result.col(result.cols-1).setTo(cv::Scalar(0));
}

inline
void
op_fliplr::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_fliplr>& in)
  {
  arma_extra_debug_sigprint();
  
  typedef typename T1::elem_type eT;
  
  const unwrap<T1>  tmp(in.m);
  const Mat<eT> X = tmp.M;
  
  if(&out != &X)
    {
    out.copy_size(X);
    
    for(uword i=0; i<X.n_cols; ++i)
      {
      out.col(i) = X.col(X.n_cols-1 - i);
      }
    }
  else
    {
    const uword N = X.n_cols / 2;
    
    for(uword i=0; i<N; ++i)
      {
      out.swap_cols(i, X.n_cols-1 - i);
      }
    }
  }

//==============================================================================
Mat
shape_model::
procrustes(const Mat &X,
           const int itol,
           const float ftol)
{
    int N = X.cols,n = X.rows/2;
    
    //remove centre of mass
    Mat P = X.clone();
    for(int i = 0; i < N; i++){
        Mat p = P.col(i);
        float mx = 0,my = 0;
        for(int j = 0; j < n; j++){mx += p.fl(2*j); my += p.fl(2*j+1);}
        mx /= n; my /= n;
        for(int j = 0; j < n; j++){p.fl(2*j) -= mx; p.fl(2*j+1) -= my;}
    }
    //optimise scale and rotation
    Mat C_old;
    for(int iter = 0; iter < itol; iter++){
        Mat C = P*Mat::ones(N,1,CV_32F)/N;
        normalize(C,C);
        if(iter > 0){if(norm(C,C_old) < ftol)break;}
        C_old = C.clone();
        for(int i = 0; i < N; i++){
            Mat R = this->rot_scale_align(P.col(i),C);
            for(int j = 0; j < n; j++){
                float x = P.fl(2*j,i),y = P.fl(2*j+1,i);
                P.fl(2*j  ,i) = R.fl(0,0)*x + R.fl(0,1)*y;
                P.fl(2*j+1,i) = R.fl(1,0)*x + R.fl(1,1)*y;
            }
        }
    }return P;
}

Mat Sharpen(Mat img)
{
    int i,j;
    Mat result;
    cvtColor(img,img,CV_BGR2GRAY);
    namedWindow("grayscale",1);
    imshow("grayscale",img);
    waitKey(2000);
    result.create(img.size(),img.type());
    for(j=1;j<img.rows-1;j++){
    uchar* previous=img.ptr<uchar>(j-1);
    uchar* current=img.ptr<uchar>(j);
    uchar* next=img.ptr<uchar>(j+1);
    uchar* output=result.ptr<uchar>(j);
    for(i=1;i<img.cols-1;i++)
    {
        *output++=cv::saturate_cast<uchar>(5*current[i]-current[i-1]-current[i+1]-previous[i]-next[i]);//saturate_cast makes resuting negative
                                                                                                  // pixel value 0 nd values above 255 255
    }
    }
    result.row(0).setTo(cv::Scalar(0));
    result.col(0).setTo(cv::Scalar(0));
    result.row(result.rows-1).setTo(cv::Scalar(0));
    result.col(result.cols-1).setTo(cv::Scalar(0));
    return result;

}

//对图像进行锐化处理
void ncu::sharpen(cv::Mat& img_in, Mat& img_out)
{
	img_out.create(img_in.size(), img_in.type());

	//处理边界内部的像素点，图像最外围的像素点应该额外处理
	for (int row = 1; row < img_in.rows - 1; row++)
	{
		//前一行像素点
		const uchar* previous = img_in.ptr<const uchar>(row - 1);
		//待处理的当前行
		const uchar* current = img_in.ptr<const uchar>(row);
		//下一行
		const uchar* next = img_in.ptr<const uchar>(row + 1);

		uchar *output = img_out.ptr<uchar>(row + 1);
		int ch = img_in.channels();
		int starts = ch;
		int ends = (img_in.cols - 1)*ch;
		for (int col = starts; col < ends; col++)
		{
			//输出图像的遍历指针与当前行的指针同步递增，以每行的每一个像素点的每一个通道值为一个递增量，
			//因为要考虑到图像的通道数
			*output++ = saturate_cast<uchar>(5 * current[col] - current[col - ch] - current[col + ch] - previous[col] - next[col]);
		}
	}//结束循环
	//处理边界，外围像素点设为 0 
	img_out.row(0).setTo(Scalar::all(0));
	img_out.row(img_out.rows - 1).setTo(Scalar::all(0));
	img_out.col(0).setTo(Scalar::all(0));
	img_out.col(img_out.cols - 1).setTo(Scalar::all(0));

}

// Draw bounding boxes on top of an image.
void	showboxes( Mat &img_color, const Mat &boxes )
{
	static const Scalar TopBoxColor = CV_RGB(255,0,0);
	static	const Scalar PartBoxColor = CV_RGB(0,0,255);

	int numfilters = int( boxes.cols / 4.f );
	for( int i=numfilters-1; i>=0; i-- ){
		Mat x1s = boxes.col( 4*i );
		Mat y1s = boxes.col( 4*i+1 );
		Mat x2s = boxes.col( 4*i+2 );
		Mat y2s = boxes.col( 4*i+3 );
		// draw each object
		for( int k=0; k<x1s.rows; k++ ){
			float	x1 = x1s.at<float>(k);
			float	y1 = y1s.at<float>(k);
			float	x2 = x2s.at<float>(k);
			float	y2 = y2s.at<float>(k);
			if( x1==0 && y1==0 && x2==0 && y2==0 )
				continue;
			Point2f		UL( x1, y1 );
			Point2f		BR( x2, y2 );
			if( i>0 )
				rectangle( img_color, UL, BR, PartBoxColor );
			else
				rectangle( img_color, UL, BR, TopBoxColor );
		}
	}

}

void Sharpen(const Mat& myImage,Mat& Result)
{
    CV_Assert(myImage.depth() == CV_8U);  // 仅接受uchar图像

    Result.create(myImage.size(),myImage.type());
    const int nChannels = myImage.channels();

    for(int j = 1 ; j < myImage.rows-1; ++j)
    {
        const uchar* previous = myImage.ptr<uchar>(j - 1);
        const uchar* current  = myImage.ptr<uchar>(j    );
        const uchar* next     = myImage.ptr<uchar>(j + 1);

        uchar* output = Result.ptr<uchar>(j);

        for(int i= nChannels;i < nChannels*(myImage.cols-1); ++i)
        {
            *output++ = saturate_cast<uchar>(5*current[i]
                         -current[i-nChannels] - current[i+nChannels] - previous[i] - next[i]);
        }
    }

    Result.row(0).setTo(Scalar(0));
    Result.row(Result.rows-1).setTo(Scalar(0));
    Result.col(0).setTo(Scalar(0));
    Result.col(Result.cols-1).setTo(Scalar(0));
}

Mat normalizePCCoeff(Mat pc, float scale, float* Cx, float* Cy, float* Cz, float* MinVal, float* MaxVal)
{
  double minVal=0, maxVal=0;

  Mat x,y,z, pcn;
  pc.col(0).copyTo(x);
  pc.col(1).copyTo(y);
  pc.col(2).copyTo(z);

  float cx = (float) cv::mean(x).val[0];
  float cy = (float) cv::mean(y).val[0];
  float cz = (float) cv::mean(z).val[0];

  cv::minMaxIdx(pc, &minVal, &maxVal);

  x=x-cx;
  y=y-cy;
  z=z-cz;
  pcn.create(pc.rows, 3, CV_32FC1);
  x.copyTo(pcn.col(0));
  y.copyTo(pcn.col(1));
  z.copyTo(pcn.col(2));

  cv::minMaxIdx(pcn, &minVal, &maxVal);
  pcn=(float)scale*(pcn)/((float)maxVal-(float)minVal);

  *MinVal=(float)minVal;
  *MaxVal=(float)maxVal;
  *Cx=(float)cx;
  *Cy=(float)cy;
  *Cz=(float)cz;

  return pcn;
}

	void CGraph::getNodes(size_t start_node, size_t num_nodes, Mat& pots) const {
		if (!num_nodes) num_nodes = getNumNodes() - start_node;

		// Assertions
		DGM_ASSERT_MSG(start_node + num_nodes <= getNumNodes(), "The given ranges exceed the number of nodes(%zu)", getNumNodes());

		if (pots.empty() || pots.cols != m_nStates || pots.rows != num_nodes)
			pots = Mat(static_cast<int>(num_nodes), m_nStates, CV_32FC1);
		
		transpose(pots, pots);

#ifdef ENABLE_PPL
		int size = pots.cols;
		int rangeSize = size / (concurrency::GetProcessorCount() * 10);
		rangeSize = MAX(1, rangeSize);
		//printf("Processors: %d\n", concurrency::GetProcessorCount());
		concurrency::parallel_for(0, size, rangeSize, [start_node, size, rangeSize, &pots, this](int i) {
			Mat pot;
			for (int j = 0; (j  < rangeSize) && (i + j < size); j++)
				getNode(start_node + i + j, lvalue_cast(pots.col(i + j)));
		});
#else
		for (int n = 0; n < pots.cols; n++)
			getNode(start_node + n, lvalue_cast(pots.col(n)));
#endif
		transpose(pots, pots);
	}

int main()
{
  BenchTimer t;
  int tries = 10;
  int rep = 400000;
  typedef Matrix3f Mat;
  typedef Vector3f Vec;
  Mat A = Mat::Random(3,3);
  A = A.adjoint() * A;

  SelfAdjointEigenSolver<Mat> eig(A);
  BENCH(t, tries, rep, eig.compute(A));
  std::cout << "Eigen:  " << t.best() << "s\n";

  Mat evecs;
  Vec evals;
  BENCH(t, tries, rep, eigen33(A,evecs,evals));
  std::cout << "Direct: " << t.best() << "s\n\n";

  std::cerr << "Eigenvalue/eigenvector diffs:\n";
  std::cerr << (evals - eig.eigenvalues()).transpose() << "\n";
  for(int k=0;k<3;++k)
    if(evecs.col(k).dot(eig.eigenvectors().col(k))<0)
      evecs.col(k) = -evecs.col(k);
  std::cerr << evecs - eig.eigenvectors() << "\n\n";
}

bool estimate_Rt_fromE(const Mat3 & K1, const Mat3 & K2,
  const Mat & x1, const Mat & x2,
  const Mat3 & E, const std::vector<size_t> & vec_inliers,
  Mat3 * R, Vec3 * t)
{
  // Accumulator to find the best solution
  std::vector<size_t> f(4, 0);

  std::vector<Mat3> Es; // Essential,
  std::vector<Mat3> Rs;  // Rotation matrix.
  std::vector<Vec3> ts;  // Translation matrix.

  Es.push_back(E);
  // Recover best rotation and translation from E.
  MotionFromEssential(E, &Rs, &ts);

  //-> Test the 4 solutions will all the point
  assert(Rs.size() == 4);
  assert(ts.size() == 4);

  Mat34 P1, P2;
  Mat3 R1 = Mat3::Identity();
  Vec3 t1 = Vec3::Zero();
  P_From_KRt(K1, R1, t1, &P1);

  for (unsigned int i = 0; i < 4; ++i)
  {
    const Mat3 &R2 = Rs[i];
    const Vec3 &t2 = ts[i];
    P_From_KRt(K2, R2, t2, &P2);
    Vec3 X;

    for (size_t k = 0; k < vec_inliers.size(); ++k)
    {
      const Vec2 & x1_ = x1.col(vec_inliers[k]),
        &x2_ = x2.col(vec_inliers[k]);
      TriangulateDLT(P1, x1_, P2, x2_, &X);
      // Test if point is front to the two cameras.
      if (Depth(R1, t1, X) > 0 && Depth(R2, t2, X) > 0)
      {
        ++f[i];
      }
    }
  }
  // Check the solution:
  const std::vector<size_t>::iterator iter = max_element(f.begin(), f.end());
  if (*iter == 0)
  {
    std::cerr << std::endl << "/!\\There is no right solution,"
      << " probably intermediate results are not correct or no points"
      << " in front of both cameras" << std::endl;
    return false;
  }
  const size_t index = std::distance(f.begin(), iter);
  (*R) = Rs[index];
  (*t) = ts[index];

  return true;
}