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

void ComputeSquaredDistanceMatrix(const vnl_matrix<double>& A,
    const vnl_matrix<double>& B, vnl_matrix<double>& D) {
  int m = A.rows();
  int n = B.rows();
  //asssert(A.cols()==B.cols());
  D.set_size(m,n);
  vnl_vector<double> v_ij;
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j < n; ++j) {
      v_ij = A.get_row(i) - B.get_row(j);
      D(i,j) = v_ij.squared_magnitude();
    }
  }
}

// todo: add one more version when the model is same as ctrl_pts
// reference:  Landmark-based Image Analysis, Karl Rohr, p195
void ComputeTPSKernel(const vnl_matrix<double>& model,
    const vnl_matrix<double>& ctrl_pts,
    vnl_matrix<double>& U, vnl_matrix<double>& K) {
  int m = model.rows();
  int n = ctrl_pts.rows();
  int d = ctrl_pts.cols();
  //asssert(model.cols()==d==(2|3));
  K.set_size(n, n);
  K.fill(0);
  U.set_size(m, n);
  U.fill(0);
  double eps = 1e-006;

  vnl_vector<double> v_ij;
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j < n; ++j) {
      v_ij = model.get_row(i) - ctrl_pts.get_row(j);
      if (d == 2) {
        double r = v_ij.squared_magnitude();
        if (r > eps) {
          U(i, j) = r * log(r) / 2;
        }
      } else if (d == 3) {
        double r = v_ij.two_norm();
        U(i, j) = -r;
      }
    }
  }
  for (int i = 0; i < n; ++i) {
    for (int j = i + 1; j < n; ++j) {
      v_ij = ctrl_pts.get_row(i) - ctrl_pts.get_row(j);
      if (d == 2) {
        double r = v_ij.squared_magnitude();
        if (r > eps) {
          K(i, j) = r * log(r) / 2;
        }
      }
      else if (d == 3) {
        double r = v_ij.two_norm();
        K(i, j) = -r;
      }
    }
  }
  for (int i = 0; i < n; ++i) {
    for (int j = 0; j < i; ++j) {
      K(i, j) = K(j, i);
    }
  }

}

vnl_matrix <double> MCLR_SM::Normalize_Feature_Matrix(vnl_matrix<double> feats)
{
	mbl_stats_nd stats;

	for(int i = 0; i<feats.rows() ; ++i)
	{
		vnl_vector<double> temp_row = feats.get_row(i);
		stats.obs(temp_row);	
	}

	std_vec = stats.sd();
	mean_vec = stats.mean();
	

//The last column is the training column 
	for(int i = 0; i<feats.columns() ; ++i)
	{
		vnl_vector<double> temp_col = feats.get_column(i);
		if(std_vec(i) > 0)
		{	
			for(int j =0; j<temp_col.size() ; ++j)
				temp_col[j] = (temp_col[j] - mean_vec(i))/std_vec(i) ;
		}
	
		feats.set_column(i,temp_col);
	}

	return feats;
}

void pick_indices(const vnl_matrix<double>&dist,
    vnl_matrix<int>&pairs, const double threshold) {
  int m = dist.rows();
  int n = dist.cols();
  vnl_vector<int> row_flag, col_flag;
  col_flag.set_size(n);  col_flag.fill(0);
  row_flag.set_size(n);  row_flag.fill(0);
  std::vector<int> row_index,col_index;
  for (int i = 0; i < m; ++i) {
    double min_dist = dist.get_row(i).min_value();
    if (min_dist < threshold) {
      for (int j = 0; j < n; ++j) {
        if (dist(i,j)==min_dist && col_flag[j] == 0){
          row_index.push_back(i);
          row_flag[i] = 1;
          col_index.push_back(j);
          col_flag[j] = 1;
        }
      }
    }
  }
  pairs.set_size(2, row_index.size());
  for (int i = 0; i<pairs.cols(); ++i){
    pairs(0,i) = row_index[i];
    pairs(1,i) = col_index[i];
  }
}

vnl_matrix <double> Normalize_Feature_Matrix(vnl_matrix<double> feats)
{
	mbl_stats_nd stats;

	for(int i = 0; i<feats.rows() ; ++i)
	{
		vnl_vector<double> temp_row = feats.get_row(i);
		stats.obs(temp_row);	
	}

	vnl_vector<double> std_vec = stats.sd();
	vnl_vector<double> mean_vec = stats.mean();
	
	for(int i = 0; i<feats.columns()-3 ; ++i)
	{
		vnl_vector<double> temp_col = feats.get_column(i);
		if(std_vec(i) > 0)
		{	
			for(int j =0; j<temp_col.size() ; ++j)
				temp_col[j] = (temp_col[j] - mean_vec(i))/std_vec(i) ;
		}
	
		feats.set_column(i,temp_col);
	}

	return feats;
}

void denormalize(vnl_matrix<double>& x, const vnl_vector<double>& centroid, const double scale) {
  int n = x.rows();
  if (n==0) return;
  int d = x.cols();
  for (int i = 0; i < n; ++i) {
    x.set_row(i, x.get_row(i) * scale + centroid);
  }
}

void normalize_same(vnl_matrix<double>& x,
    vnl_vector<double>& centroid, double& scale) {
  int n = x.rows();
  if (n==0) return;
  int d = x.cols();
  for (int i = 0; i < n; ++i) {
    x.set_row(i, (x.get_row(i) - centroid) / scale);
  }
}

void f(const vnl_matrix<double>& model,
    const vnl_matrix<double>& scene, double threshold,
    vnl_matrix<double>& extracted_model,
    vnl_matrix<double>& extracted_scene) {
  vnl_matrix<double> dist;
  vnl_matrix<int> pairs;
  ComputeSquaredDistanceMatrix(model, scene, dist);
  pick_indices(dist, pairs, threshold*threshold);
  std::cout << "distance threshold : " << threshold << std::endl;
  int j, n = pairs.cols();
  int d = model.cols();
  extracted_model.set_size(n,d);
  extracted_scene.set_size(n,d);
  std::cout << "# of matched point pairs : " << n << std::endl;
  for (j=0; j<n; ++j) {
    extracted_model.set_row(j,model.get_row(pairs(0,j)));
  }
  for (j=0; j<n; ++j) {
    extracted_scene.set_row(j,scene.get_row(pairs(1,j)));
  }
}

void ComputeTPSKernelU(const vnl_matrix<double>& model,
    const vnl_matrix<double>& ctrl_pts,
    vnl_matrix<double>& U) {
  int m = model.rows();
  int n = ctrl_pts.rows();
  int d = ctrl_pts.cols();
  //asssert(model.cols()==d==(2|3));
  //K.set_size(n, n);
  //K.fill(0);
  U.set_size(m, n);
  U.fill(0);
  double eps = 1e-006;

  vnl_vector<double> v_ij;
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j < n; ++j) 
	{
      v_ij = model.get_row(i) - ctrl_pts.get_row(j);
      double r = v_ij.two_norm();
      U(i, j) = -r;
    }
  }
}

void ExtractMatchingPairs(
    const vnl_matrix<T>& model,
    const vnl_matrix<T>& scene,
    const T& threshold,
    vnl_matrix<T>& extracted_model,
    vnl_matrix<T>& extracted_scene) {
  vnl_matrix<T> dist;
  vnl_matrix<int> pairs;
  ComputeSquaredDistanceMatrix<T>(model, scene, dist);
  PickIndices<T>(dist, pairs, threshold*threshold);
  std::cout << "distance threshold : " << threshold << std::endl;
  int n = pairs.cols();
  int d = model.cols();
  extracted_model.set_size(n, d);
  extracted_scene.set_size(n, d);
  std::cout << "# of matched point pairs : " << n << std::endl;
  for (int j = 0; j < n; ++j) {
    extracted_model.set_row(j,model.get_row(pairs(0, j)));
  }
  for (int j = 0; j < n; ++j) {
    extracted_scene.set_row(j,scene.get_row(pairs(1, j)));
  }
}

void itk::BiExpFitFunctor::operator()(vnl_matrix<double> & newSignal,const vnl_matrix<double> & SignalMatrix, const double & S0)
{

  vnl_vector<double> initalGuess(3);
  // initialize Least Squres Function
  // SignalMatrix.cols() defines the number of shells points
  lestSquaresFunction model(SignalMatrix.cols());
  model.set_bvalues(m_BValueList);// set BValue Vector e.g.: [1000, 2000, 3000] <- shell b Values

  // initialize Levenberg Marquardt
  vnl_levenberg_marquardt minimizer(model);
  minimizer.set_max_function_evals(1000);   // Iterations
  minimizer.set_f_tolerance(1e-10);        // Function tolerance

  // for each Direction calculate LSF Coeffs ADC & AKC
  for(unsigned int i = 0 ; i < SignalMatrix.rows(); i++)
  {
    model.set_measurements(SignalMatrix.get_row(i));
    model.set_reference_measurement(S0);

    initalGuess.put(0, 0.f); // ADC_slow
    initalGuess.put(1, 0.009f); // ADC_fast
    initalGuess.put(2, 0.7f); // lambda

    // start Levenberg-Marquardt
    minimizer.minimize_without_gradient(initalGuess);

    const double & ADC_slow = initalGuess.get(0);
    const double & ADC_fast = initalGuess.get(1);
    const double & lambda = initalGuess(2);

    newSignal.put(i, 0, S0 * (lambda * std::exp(-m_TargetBvalue * ADC_slow) + (1-lambda)* std::exp(-m_TargetBvalue * ADC_fast)));
    newSignal.put(i, 1, minimizer.get_end_error()); // RMS Error

    //OUTPUT FOR EVALUATION
    /*std::cout << std::scientific << std::setprecision(5)
              << ADC_slow   << ","                        // lambda
              << ADC_fast   << ","                        // alpha
              << lambda     << ","                        // lambda
              << S0         << ","                        // S0 value
              << minimizer.get_end_error() << ",";      // End error
    for(unsigned int j = 0; j < SignalMatrix.get_row(i).size(); j++ ){
      std::cout << std::scientific << std::setprecision(5) << SignalMatrix.get_row(i)[j];    // S_n Values corresponding to shell 1 to shell n
      if(j != SignalMatrix.get_row(i).size()-1) std::cout << ",";
    }
    std::cout << std::endl;*/
  }

}

void normalize(vnl_matrix<double>& x,
    vnl_vector<double>& centroid, double& scale) {
  int n = x.rows();
  if (n==0) return;
  int d = x.cols();
  centroid.set_size(d);

  vnl_vector<double> col;
  for (int i = 0; i < d; ++i) {
    col = x.get_column(i);
    centroid(i) = col.mean();
  }
  for (int i = 0; i < n; ++i) {
    x.set_row(i, x.get_row(i) - centroid);
  }
  scale = x.frobenius_norm() / sqrt(double(n));
  x = x / scale;
}

// x contains the training features
// y contains the training labels
void MCLR_SM::Initialize(vnl_matrix<double> data,double c,vnl_vector<double> classes, std::string str,vtkSmartPointer<vtkTable> table )
{
	int train_counter =0;
	int test_counter = 0;
	validation = str;
	
	// test_table contains the table after the queried samples are removed from the original table
	// It is updated in the Update_Train_data 
	test_table = table;
	
	for(int i = 0; i< classes.size() ; ++i)
	{
		if(classes(i) == -1)
			test_counter++;
	}
	
	train_data.set_size(data.rows()-test_counter,data.cols());
	test_data.set_size(test_counter,data.cols());
	
	test_counter = 0;

	for(int i= 0;i<data.rows();++i)
	{
		if(classes(i)!=-1)
		{	
			train_data.set_row(train_counter,data.get_row(i));
			++train_counter;
			test_table->RemoveRow(i);
		}
		else
		{
			test_data.set_row(test_counter,data.get_row(i));
			++test_counter;
		}
	}
	

	y.set_size(train_data.rows());

	// Transpose : Features->Rows ; Samples -> Columns
	x = train_data.transpose();// Feature matrix of training samples only ; does not contain unlabeled sample data
	
	int counter = 0;
	for(int i = 0; i < classes.size(); ++i)
	{
		if(classes(i)!=-1)
		{
			y(counter) = classes(i);// class value ; 
			counter++;
		}
	}
	


	//m.method = "direct"; // No kernel is the default
	m.sparsity_control = c; // greater the value , more is the sparsity
	m.kstd_level = 1; // kernel sigma val

	// Stop Condition
	// Used in active label selection
	stop_cond.set_size(2);
	stop_cond[0] = 1e-5;
	stop_cond[1] = 0.01;
	delta = 1e-9; // used for approximating |w| 
	
	no_of_features = x.rows();
	// We assume that the classes are numbered from 1 to "n" for n classes
	// so y.maxval - .minval +1 gives number of classes
	no_of_classes = y.max_value()-y.min_value()+1;
	class_vector.set_size(no_of_classes);

	m.w.set_size(no_of_features+1,no_of_classes);
	m.w.fill(0);

	gradient_w.set_size(no_of_features+1,no_of_classes);
	hessian.set_size((no_of_features+1)*(no_of_classes),(no_of_features+1)*(no_of_classes)); /// Check the size again!!!
	
	// g :the objective function value
	g = 0;
	// Class vector
	for(int i=1;i<=no_of_classes;++i)
	{
		class_vector.put(i-1,i); 
	}
	
	diff_info_3_it.set_size(3);
	info_3_it.set_size(3);

}