C++ vnl_matrix::empty方法代码示例

void compute_P(const vnl_matrix<double>& x,const vnl_matrix<double>& y, vnl_matrix<double>& P, double &E, double sigma, int outliers) {
  double k;
  k = -2*sigma*sigma;

  //P.set_size(m,n); P.fill(0);
  //vnl_vector<double> v_ij;

  vnl_vector<double> column_sum;
  int m = x.rows();
  int s = y.rows();
  int d = x.cols();
  column_sum.set_size(s);
  column_sum.fill(0);
  double outlier_term = outliers*pow((2*sigma*sigma*3.1415926),0.5*d);
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j < s; ++j) {
      double r = 0;
      for (int t = 0; t < d; ++t) {
        r += (x(i,t) - y(j,t))*(x(i,t) - y(j,t));
      }
      P(i,j) = exp(r/k);
      column_sum[j]+=P(i,j);
    }
  }


  if (outliers!=0) {
    for (int i = 0; i < s; ++i)
      column_sum[i] += outlier_term;
  }
  if (column_sum.min_value()>(1e-12)) {
    E = 0;
    for (int i = 0; i < s; ++i) {
      for (int j = 0; j < m; ++j){
        P(j,i) = P(j,i)/column_sum[i];
      }
      E -= log(column_sum[i]);
    }
    //vcl_cerr< < s;
    //vcl_cerr<<P.get_column(10);
  }
  else {
    P.empty();
  }
}

void RegPowell::computeIterativeRegistration( int nmax,double epsit,MRI * mriS, MRI* mriT, const vnl_matrix < double > & m , double iscaleinit)
// retruns 4x4 matrix and iscale value
{
  if (!mriS) mriS = mri_source;
  if (!mriT) mriT = mri_target;

  assert (mriS && mriT);

  tocurrent = this; // so that we can access this from static cost function

  pair < vnl_matrix_fixed < double, 4, 4> , double > fmd(vnl_matrix_fixed < double, 4 , 4> () ,iscaleinit);

  // check if mi (inital transform) is passed
  if (!m.empty()) fmd.first = m;
  else if (!Minit.empty()) fmd.first = Minit;
  else fmd.first = initializeTransform(mriS,mriT) ;

  if (debug > 0)
  {
    cout << "   - initial transform:\n" ;
    //MatrixPrintFmt(stdout,"% 2.8f",fmd.first);
		cout << fmd.first << endl;
  }


  // here maybe better to symmetrically warp both images SQRT(M)
  // this keeps the problem symmetric
  cout << "   - warping source and target (sqrt)" << endl;
  //if (mh1) MatrixFree(&mh1);
  mh1 = MyMatrix::MatrixSqrt(fmd.first);
  // do not just assume m = mh*mh, rather m = mh2 * mh
  // for transforming target we need mh2^-1 = mh * m^-1
  //MATRIX * mi  = MatrixInverse(fmd.first,NULL);
  //MATRIX * mhi = MatrixMultiply(mh1,mi,NULL);
  vnl_matrix_fixed < double, 4, 4 > mhi = mh1 * vnl_inverse(fmd.first);
	
  //set static
  //if (mh2) MatrixFree(&mh2);
  //mh2 = MatrixInverse(mhi,NULL); // M = mh2 * mh1
  mh2 = vnl_inverse(mhi); // M = mh2 * mh1

  //if (mri_Swarp) MRIfree(&mri_Swarp);
  //mri_Swarp = MRIclone(mriS,NULL);
  //mri_Swarp = MRIlinearTransform(mriS,mri_Swarp, mh);
  //if (mri_Twarp) MRIfree(&mri_Twarp);
  //mri_Twarp = MRIclone(mriS,NULL); // bring them to same space (just use src geometry)
  //mri_Twarp = MRIlinearTransform(mriT,mri_Twarp, mhi);

 // //MatrixFree(&mh);
  //MatrixFree(&mhi);
 // MatrixFree(&mi);

  // adjust intensity later in first powell call

  scf = mriS;
  tcf = mriT;

  // create parameter vector:
  pcount = 3; // transolny
  if (rigid) pcount = 6;
  else pcount = 12;
  if (pcount==3) assert(transonly);
  if (iscale) pcount++;


  // compute Registration
  cout << "   - compute new registration ( " << pcount << " params )" << endl;
  float fret, fstart, min_sse;
  float* p   = ::vector(1, pcount+1) ;
  float** xi = ::matrix(1, pcount+1, 1, pcount+1) ;
  float tol = 1e-5; //-8
  int maxiter = 36;
  int iter;

  for (int i = 1;i<=pcount;i++)
  {
    p[i] = 0.0;

    for (int j = i;j<=pcount;j++)
    {
      xi[i][j] = 0.0;
      xi[j][i] = 0.0;
    }
    xi[i][i] = 1.0;
  }
  if (iscale) p[pcount] = iscaleinit;

  min_sse = costFunction(p) ;
  cout << " min_sse: " << min_sse << endl;
  icount = 0;
  //OpenPowell(p, xi, pcount, tol, &iter, &fret, costFunction);
  OpenPowell2(p, xi, pcount, tol,tol,maxiter, &iter, &fret, costFunction);
  cout << endl << "best alignment initial powell: " << fret << " (" << iter << " steps)" << endl;

  int count = 0;
  do
  {
    count++;
    // reinitialize powell directions
    for (int r = 1 ; r <= pcount ; r++)
//.........这里部分代码省略.........