C++ Matrix::row方法代码示例

    void solve_nonlinear () {
      for (size_t i = 0; i < signals.rows(); ++i) {
        const Math::Vector<cost_value_type> signal (signals.row(i));
        Math::Vector<cost_value_type> values (tensors.row(i));

        cost.set_voxel (&signal, &values);

        Math::Vector<cost_value_type> x (cost.size());
        cost.init (x);
        //Math::check_function_gradient (cost, x, 1e-10, true);

        Math::GradientDescent<Cost> optim (cost);
        try { optim.run (10000, 1e-8); }
        catch (Exception& E) {
          E.display();
        }

        //x = optim.state();
        //Math::check_function_gradient (cost, x, 1e-10, true);

        cost.get_values (values, optim.state());
      }
    }

void run () 
{
  try {
    Math::Matrix<value_type> directions = DWI::Directions::load_cartesian<value_type> (argument[0]);
    report (str(argument[0]), directions);
  }
  catch (Exception& E) {
    Math::Matrix<value_type> directions (str(argument[0]));
    DWI::normalise_grad (directions);
    if (directions.columns() < 3) 
      throw Exception ("unexpected matrix size for DW scheme \"" + str(argument[0]) + "\"");

    print (str(argument[0]) + " [ " + str(directions.rows()) + " volumes ]\n");
    DWI::Shells shells (directions);

    for (size_t n = 0; n < shells.count(); ++n) {
      Math::Matrix<value_type> subset (shells[n].count(), 3);
      for (size_t i = 0; i < subset.rows(); ++i)
        subset.row(i) = directions.row(shells[n].get_volumes()[i]).sub(0,3);
      report ("\nb = " + str(shells[n].get_mean()), subset);
    }
  }
}

void run () 
{
  Math::Matrix<value_type> directions = DWI::Directions::load_cartesian<value_type> (argument[0]);

  size_t num_permutations = 1e8;
  Options opt = get_options ("permutations");
  if (opt.size())
    num_permutations = opt[0][0];

  Shared eddy_shared (directions, num_permutations);
  Thread::run (Thread::multi (Processor (eddy_shared)), "eval thread");

  auto& signs = eddy_shared.get_best_signs();

  for (size_t n = 0; n < directions.rows(); ++n) 
    if (signs[n] < 0)
      directions.row(n) *= -1.0;

  bool cartesian = get_options("cartesian").size();
  DWI::Directions::save (directions, argument[1], cartesian);
}

void report (const std::string& title, const Math::Matrix<value_type>& directions) 
{
  std::vector<value_type> NN_bipolar (directions.rows(), 0.0);
  std::vector<value_type> NN_unipolar (directions.rows(), 0.0);

  std::vector<value_type> E_bipolar (directions.rows(), 0.0);
  std::vector<value_type> E_unipolar (directions.rows(), 0.0);

  for (size_t i = 0; i < directions.rows()-1; ++i) {
    for (size_t j = i+1; j < directions.rows(); ++j) {
      value_type cos_angle = Math::dot (directions.row(i).sub(0,3), directions.row(j).sub(0,3));
      NN_unipolar[i] = std::max (NN_unipolar[i], cos_angle);
      NN_unipolar[j] = std::max (NN_unipolar[j], cos_angle);
      cos_angle = std::abs(cos_angle);
      NN_bipolar[i] = std::max (NN_bipolar[i], cos_angle);
      NN_bipolar[j] = std::max (NN_bipolar[j], cos_angle);

      value_type E = 
        Math::pow2 (directions(i,0) - directions(j,0)) + 
        Math::pow2 (directions(i,1) - directions(j,1)) + 
        Math::pow2 (directions(i,2) - directions(j,2));
      E = value_type (1.0) / E;

      E_unipolar[i] += E;
      E_unipolar[j] += E;

      value_type E2 = 
        Math::pow2 (directions(i,0) + directions(j,0)) + 
        Math::pow2 (directions(i,1) + directions(j,1)) + 
        Math::pow2 (directions(i,2) + directions(j,2));
      E += value_type (1.0) / E2;

      E_bipolar[i] += E;
      E_bipolar[j] += E;

    }
  }




  auto report_NN = [](const std::vector<value_type>& NN) {
    value_type min = std::numeric_limits<value_type>::max();
    value_type mean = 0.0;
    value_type max = 0.0;
    for (auto a : NN) {
      a = (180.0/Math::pi) * std::acos (a);
      mean += a;
      min = std::min (min, a);
      max = std::max (max, a);
    }
    mean /= NN.size();

    print ("    nearest-neighbour angles: mean = " + str(mean) + ", range [ " + str(min) + " - " + str(max) + " ]\n");
  };



  auto report_E = [](const std::vector<value_type>& E) {
    value_type min = std::numeric_limits<value_type>::max();
    value_type total = 0.0;
    value_type max = 0.0;
    for (auto e : E) {
      total += e;
      min = std::min (min, e);
      max = std::max (max, e);
    }
    print ("    energy: total = " + str(total) + ", mean = " + str(total/E.size()) + ", range [ " + str(min) + " - " + str(max) + " ]\n");
  };




  print (title + " [ " + str(directions.rows()) + " directions ]\n\n");

  print ("  Bipolar electrostatic repulsion model:\n");
  report_NN (NN_bipolar);
  report_E (E_bipolar);

  print ("\n  Unipolar electrostatic repulsion model:\n");
  report_NN (NN_unipolar);
  report_E (E_unipolar);

  std::string lmax_results;
  for (size_t lmax = 2; lmax <= Math::SH::LforN (directions.rows()); lmax += 2) 
    lmax_results += " " + str(DWI::condition_number_for_lmax (directions, lmax));
  print ("\n  Spherical Harmonic fit:\n    condition numbers for lmax = " + str(2) + " -> " 
      + str(Math::SH::LforN (directions.rows())) + ":" + lmax_results + "\n\n");
}

void run() {
  Image::BufferPreload<float> data_in (argument[0], Image::Stride::contiguous_along_axis (3));
  auto voxel_in = data_in.voxel();

  Math::Matrix<value_type> grad (DWI::get_valid_DW_scheme<float> (data_in));

  // Want to support non-shell-like data if it's just a straight extraction
  //   of all dwis or all bzeros i.e. don't initialise the Shells class
  std::vector<size_t> volumes;
  bool bzero = get_options ("bzero").size();
  Options opt = get_options ("shell");
  if (opt.size()) {
    DWI::Shells shells (grad);
    shells.select_shells (false, false);
    for (size_t s = 0; s != shells.count(); ++s) {
      DEBUG ("Including data from shell b=" + str(shells[s].get_mean()) + " +- " + str(shells[s].get_stdev()));
      for (std::vector<size_t>::const_iterator v = shells[s].get_volumes().begin(); v != shells[s].get_volumes().end(); ++v)
        volumes.push_back (*v);
    }
    // Remove DW information from header if b=0 is the only 'shell' selected
    bzero = (shells.count() == 1 && shells[0].is_bzero());
  } else {
    const float bzero_threshold = File::Config::get_float ("BValueThreshold", 10.0);
    for (size_t row = 0; row != grad.rows(); ++row) {
      if ((bzero && (grad (row, 3) < bzero_threshold)) || (!bzero && (grad (row, 3) > bzero_threshold)))
        volumes.push_back (row);
    }
  }

  if (volumes.empty())
    throw Exception ("No " + str(bzero ? "b=0" : "dwi") + " volumes present");

  std::sort (volumes.begin(), volumes.end());

  Image::Header header (data_in);

  if (volumes.size() == 1)
    header.set_ndim (3);
  else
    header.dim (3) = volumes.size();

  Math::Matrix<value_type> new_grad (volumes.size(), grad.columns());
  for (size_t i = 0; i < volumes.size(); i++)
    new_grad.row (i) = grad.row (volumes[i]);
  header.DW_scheme() = new_grad;

  Image::Buffer<value_type> data_out (argument[1], header);
  auto voxel_out = data_out.voxel();

  Image::Loop outer ("extracting volumes...", 0, 3);

  if (voxel_out.ndim() == 4) {

    for (auto i = outer (voxel_out, voxel_in); i; ++i) {
      for (size_t i = 0; i < volumes.size(); i++) {
        voxel_in[3] = volumes[i];
        voxel_out[3] = i;
        voxel_out.value() = voxel_in.value();
      }
    }

  } else {

    const size_t volume = volumes[0];
    for (auto i = outer (voxel_out, voxel_in); i; ++i) {
      voxel_in[3] = volume;
      voxel_out.value() = voxel_in.value();
    }

  }


}