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

	DenseKernel(const MatrixXf & f, KernelType ktype, NormalizationType ntype):f_(f), ktype_(ktype), ntype_(ntype) {
		if (ktype_ == DIAG_KERNEL)
			parameters_ = VectorXf::Ones( f.rows() );
		else if( ktype == FULL_KERNEL )
			parameters_ = MatrixXf::Identity( f.rows(), f.rows() );
		initLattice( f );
	}

	MatrixXf featureGradient( const MatrixXf & a, const MatrixXf & b ) const {
		if (ntype_ == NO_NORMALIZATION )
			return kernelGradient( a, b );
		else if (ntype_ == NORMALIZE_SYMMETRIC ) {
			MatrixXf fa = lattice_.compute( a*norm_.asDiagonal(), true );
			MatrixXf fb = lattice_.compute( b*norm_.asDiagonal() );
			MatrixXf ones = MatrixXf::Ones( a.rows(), a.cols() );
			VectorXf norm3 = norm_.array()*norm_.array()*norm_.array();
			MatrixXf r = kernelGradient( 0.5*( a.array()*fb.array() + fa.array()*b.array() ).matrix()*norm3.asDiagonal(), ones );
			return - r + kernelGradient( a*norm_.asDiagonal(), b*norm_.asDiagonal() );
		}
		else if (ntype_ == NORMALIZE_AFTER ) {
			MatrixXf fb = lattice_.compute( b );
			
			MatrixXf ones = MatrixXf::Ones( a.rows(), a.cols() );
			VectorXf norm2 = norm_.array()*norm_.array();
			MatrixXf r = kernelGradient( ( a.array()*fb.array() ).matrix()*norm2.asDiagonal(), ones );
			return - r + kernelGradient( a*norm_.asDiagonal(), b );
		}
		else /*if (ntype_ == NORMALIZE_BEFORE )*/ {
			MatrixXf fa = lattice_.compute( a, true );
			
			MatrixXf ones = MatrixXf::Ones( a.rows(), a.cols() );
			VectorXf norm2 = norm_.array()*norm_.array();
			MatrixXf r = kernelGradient( ( fa.array()*b.array() ).matrix()*norm2.asDiagonal(), ones );
			return -r+kernelGradient( a, b*norm_.asDiagonal() );
		}
	}

// normalizeMatch respect to "In defense of eight point algorithm"
void normalizeMatch(MatrixXf &mat, Matrix3f &T1, Matrix3f &T2) {
  MatrixXf pts1 = mat.leftCols<3>();
  MatrixXf pts2 = mat.block(0, 3, mat.rows(), 3);
  normalizePts(pts1, T1);
  normalizePts(pts2, T2);
  mat.leftCols<3>() = pts1;
  mat.block(0, 3, mat.rows(), 3) = pts2;
}

void Permutohedral::compute ( MatrixXf & out, const MatrixXf & in, bool reverse ) const
{
    if( out.cols() != in.cols() || out.rows() != in.rows() )
        out = 0*in;
    if( in.rows() <= 2 )
        seqCompute( out.data(), in.data(), in.rows(), reverse );
    else
        sseCompute( out.data(), in.data(), in.rows(), reverse );
}

void blas_gemm(const MatrixXf& a, const MatrixXf& b, MatrixXf& c)
{
  int M = c.rows(); int N = c.cols(); int K = a.cols();
  int lda = a.rows(); int ldb = b.rows(); int ldc = c.rows();

  sgemm_(&notrans,&notrans,&M,&N,&K,&fone,
         const_cast<float*>(a.data()),&lda,
         const_cast<float*>(b.data()),&ldb,&fone,
         c.data(),&ldc);
}

QPainterPath Layouter::mat2Path( const MatrixXf& pntMat )
{
	QPainterPath path;
	if (pntMat.rows() <= 0 || pntMat.cols() != 2)
		return path;
	path.moveTo(pntMat(0,0), pntMat(0,1));
	for (int i = 1; i < pntMat.rows(); ++i)
		path.lineTo(pntMat(i,0), pntMat(i,1));
	return path;
}

bool singleModelRANSAC(const MatrixXf &data, int M, MatrixXf &inlier) {
  int maxdegen = 10;
  int dataSize = data.rows();
  int psize = 4;
  MatrixXf x1 = data.block(0, 0, data.rows(), 3);
  MatrixXf x2 = data.block(0, 3, data.rows(), 3);
  vector<int> sample;
  MatrixXf pts1(4, 3);
  MatrixXf pts2(4, 3);
  int maxInlier = -1;
  MatrixXf bestResidue;
  for (int m = 0; m < M; m++) {
    int degencount = 0;
    int isdegen = 1;
    while (isdegen==1 && degencount < maxdegen) {
      degencount ++;
      RandomSampling(psize, dataSize, sample);
      for (int i = 0; i < psize; i++) {
        pts1.row(i) = x1.row(sample[i]);
        pts2.row(i) = x2.row(sample[i]);
      }
      if (sampleValidTest(pts1, pts2))
          isdegen = 0;
    }
    if (isdegen) {
      cout << "Cannot find valid p-subset" << endl;
      return false;
    }
    Matrix3f local_H;
    MatrixXf local_A;
    fitHomography(pts1, pts2, local_H, local_A);

    MatrixXf residue;
    computeHomographyResidue(x1, x2, local_H, residue);
    int inlierCount = (residue.array() < THRESHOLD).count();
    if (inlierCount > maxInlier) {
      maxInlier = inlierCount;
      bestResidue = residue;
    }
  }
  inlier.resize(maxInlier, data.cols());
  int transferCounter = 0;
  for (int i = 0; i < dataSize; i++) {
    if (bestResidue(i) < THRESHOLD) {
      inlier.row(transferCounter) = data.row(i);
      transferCounter++;
    }
  }
  if (transferCounter != maxInlier) {
    cout << "RANSAC result size does not match!!!!" << endl;
    return false;
  }
  return true;
}

void SortPoints(MatrixXf &x2d)
{
  // the array is formed by [x1 y1 1,
//                             x2 y2  ]
  
  
  MatrixXf x2d_aux=MatrixXf::Zero(x2d.rows(),3);

  //////////////////////////////////////////////////////////////
  //first sort point in y
  ///////////////////////////////////////////////////////////////
  
  int i,j;
  for (i =0; i<x2d.rows();i++){
    for(j = i+1; j < x2d.rows(); j ++) {
        if(x2d(j,1) < x2d(i,1)) {
               float temp_x = x2d(i,0);
	       float temp_y = x2d(i,1);
               x2d(0,i) = x2d(0,j);
	       x2d(1,i) = x2d(1,j);
	       
	       x2d(0,j) = temp_x;
	       x2d(1,j) = temp_y;
               
               
	}
    
    }
  }
 
 /////////////////////////////////////////
  //now order in X
 
  for (i =0; i<x2d.rows();i++){
    for(j = i+1; j < x2d.rows(); j ++) {
        if(x2d(j,0) < x2d(i,0)) {
               float temp_x = x2d(i,0);
	       float temp_y = x2d(i,1);
               x2d(0,i) = x2d(0,j);
	       x2d(1,i) = x2d(1,j);
	       
	       x2d(0,j) = temp_x;
	       x2d(1,j) = temp_y;
               
               
	}
    
    }
  }
  
  
}

vector<vector<float> > applyPCAtoVector2D(vector<vector<float> > &descriptorValues, MatrixXf &eigen_vects)
{
    MatrixXf datapoints(descriptorValues.size(),descriptorValues[0].size());
    for (int i = 0; i < descriptorValues.size(); ++i)
        for (int j = 0; j < descriptorValues[0].size(); ++j)
            datapoints(i, j) = descriptorValues[i][j];
    MatrixXf reduceddatapnts = pca::transformPointMatrix(datapoints, eigen_vects);
    vector<vector<float> > retfeatvects(reduceddatapnts.rows(), vector<float>(reduceddatapnts.cols()));
    for (int i = 0; i < reduceddatapnts.rows(); ++i)
        for (int j = 0; j < reduceddatapnts.cols(); ++j)
            retfeatvects[i][j] = reduceddatapnts(i,j);
    return retfeatvects;
}

void Neuromag::run()
{
    MatrixXf matValue;

    qint32 size = 0;

    while(m_bIsRunning) {
        if(m_pRawMatrixBuffer_In) {
            //pop matrix
            matValue = m_pRawMatrixBuffer_In->pop();

            //Write raw data to fif file
            if(m_bWriteToFile) {
                size += matValue.rows()*matValue.cols() * 4;

                if(size > MAX_DATA_LEN) {
                    size = 0;
                    this->splitRecordingFile();
                }

                m_mutex.lock();
                if(m_pOutfid) {
                    m_pOutfid->write_raw_buffer(matValue.cast<double>());
                }
                m_mutex.unlock();
            } else {
                size = 0;
            }

            if(m_pRTMSA_Neuromag) {
                m_pRTMSA_Neuromag->data()->setValue(this->calibrate(matValue));
            }
        }
    }
}

VectorXf EMclustering::logsumexp(MatrixXf x, int dim)
{
	int r = x.rows();
	int c = x.cols();

	VectorXf y(r);
	MatrixXf tmp1(r,c);
	VectorXf tmp2(r);
	VectorXf s(r);

	y = x.rowwise().maxCoeff();//cerr<<"y"<<y<<endl<<endl;
	x = x.colwise() - y;	
	//cerr<<"x"<<x<<endl<<endl;
	tmp1 = x.array().exp();	
	//cerr<<"t"<<tmp1<<endl<<endl;
	tmp2 = tmp1.rowwise().sum();	
	//cerr<<"t"<<tmp2<<endl<<endl;
	s = y.array() + tmp2.array().log();

	for(int i=0;i<s.size();i++)
	{
		if(!isfinite(s(i)))
		{
			s(i) = y(i);
		}
	}

	y.resize(0);
	tmp1.resize(0,0);
	tmp2.resize(0);
	
	return s;
}

VectorXf param_sensitivity_widget::computeSensitivity(
    MatrixXf &parameterMatrix, VectorXf &responseVector)
{
    MatrixXf Ctemp = parameterMatrix.transpose()*parameterMatrix;
    MatrixXf C;
    C = Ctemp.inverse();

    VectorXf b = C*parameterMatrix.transpose()*responseVector;

    VectorXf Y_hat = parameterMatrix*b;

    int p = b.rows();

    VectorXf sigma2Vec = responseVector-Y_hat;

    float sigma2 = sigma2Vec.squaredNorm();
    sigma2= sigma2/(parameterMatrix.rows() - p);

    Ctemp = C*sigma2;

    MatrixXf denominator = Ctemp.diagonal();

    // Do element-wise division
    VectorXf t = b;
    for (int i = 0; i < b.rows(); i++)
    {
        t(i) = abs(b(i)/sqrt(denominator(i)));
    }

    return t;
}

void FiffSimulator::doContinousHPI(MatrixXf& matData)
{
    //This only works with babyMEG HPI channels 400 ... 407
    if(m_pFiffInfo && m_pHPIWidget && matData.rows() >= 407) {
        if(m_pHPIWidget->wasLastFitOk()) {
            // Load device to head transformation matrix from Fiff info
            QMatrix3x3 rot;

            for(int ir = 0; ir < 3; ir++) {
                for(int ic = 0; ic < 3; ic++) {
                    rot(ir,ic) = m_pFiffInfo->dev_head_t.trans(ir,ic);
                }
            }

            QQuaternion quatHPI = QQuaternion::fromRotationMatrix(rot);

            // Write rotation quaternion to HPI Ch #1~3
            matData.row(401) = MatrixXf::Constant(1,matData.cols(), quatHPI.x());
            matData.row(402) = MatrixXf::Constant(1,matData.cols(), quatHPI.y());
            matData.row(403) = MatrixXf::Constant(1,matData.cols(), quatHPI.z());

            // Write translation vector to HPI Ch #4~6
            matData.row(404) = MatrixXf::Constant(1,matData.cols(), m_pFiffInfo->dev_head_t.trans(0,3));
            matData.row(405) = MatrixXf::Constant(1,matData.cols(), m_pFiffInfo->dev_head_t.trans(1,3));
            matData.row(406) = MatrixXf::Constant(1,matData.cols(), m_pFiffInfo->dev_head_t.trans(2,3));

            // Write GOF to HPI Ch #7
            // Write goodness of fit (GOF)to HPI Ch #7
            float dpfitError = 0.0;
            float GOF = 1 - dpfitError;
            matData.row(407) = MatrixXf::Constant(1,matData.cols(), GOF);
        }
    }
}

double IntersectionOverUnion::evaluate( MatrixXf & d_mul_Q, const MatrixXf & Q ) const {
    assert( gt_.rows() == Q.cols() );
    const int N = Q.cols(), M = Q.rows();
    d_mul_Q = 0*Q;

    VectorXd in(M), un(M);
    in.fill(0.f);
    un.fill(1e-20);
    for( int i=0; i<N; i++ ) {
        if( 0 <= gt_[i] && gt_[i] < M ) {
            in[ gt_[i] ] += Q(gt_[i],i);
            un[ gt_[i] ] += 1;
            for( int l=0; l<M; l++ )
                if( l!=gt_[i] )
                    un[ l ] += Q(l,i);
        }
    }
    for( int i=0; i<N; i++ )
        if( 0 <= gt_[i] && gt_[i] < M ) {
            for( int l=0; l<M; l++ )
                if( l==gt_[i] )
                    d_mul_Q(l,i) = Q(l,i) / (un[l]*M);
                else
                    d_mul_Q(l,i) = - Q(l,i) * in[l] / ( un[l] * un[l] * M);
        }
    return (in.array()/un.array()).sum()/M;
}

float get_te_cost(int row, int col, int i, const MatrixXf & cost, const ExternNDArrayf & gradient_img) {
    if ((row + i < 0) || (row + i >= cost.rows())) {
        return INFINITY;
    } else {
        return (col == 0 ? 0 : cost(row+i, col-1)) + gradient_img(row, col);
    }
}