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

 DWI2QBI (const Math::Matrix<value_type>& FRT_SHT, Math::Matrix<value_type>& normalise_SHT, const DWI::Shells& shells) :
   FRT_SHT (FRT_SHT), 
   normalise_SHT (normalise_SHT), 
   shells (shells),
   dwi (FRT_SHT.columns()),
   qbi (FRT_SHT.rows()),
   amps (normalise ? normalise_SHT.rows() : 0) { }

void run ()
{

  DWI::Tractography::Properties properties;

  DWI::Tractography::Writer<> writer (argument.back(), properties);

  for (size_t n = 0; n < argument.size()-1; n++) {
    Math::Matrix<float> M;
    try {
      M.load (argument[n]);
      if (M.columns() != 3)
        throw Exception ("file \"" + argument[n] + "\" does not contain 3 columns - ignored");

      DWI::Tractography::Streamline<float> tck (M.rows());
      for (size_t i = 0; i < M.rows(); i++) {
        tck[i].set (M (i,0), M (i,1), M (i,2));
      }
      writer (tck);
      writer.total_count++;
    }
    catch (Exception) { }
  }

}

void
AAKR::add(Math::Matrix v)
{
    if (dataSize() == 0)
        throw std::runtime_error("unable to add: data window size is undefined.");

    if (v.rows() != 1)
        throw std::runtime_error("unable to add: new sample is not a row vector.");

    if (sampleSize() == 0)
        m_data.resize(dataSize(), v.columns());

    if ((unsigned)v.columns() != sampleSize())
        throw std::runtime_error("unable to add: sample size does not match.");

    // Write to the data set.
    m_data.set(m_index, m_index, 0, sampleSize() - 1, v);

    // Increment data set index.
    increment();
}

void run()
{
  InputBufferType dwi_buffer (argument[0], Image::Stride::contiguous_along_axis (3));
  Math::Matrix<cost_value_type> grad = DWI::get_valid_DW_scheme<cost_value_type> (dwi_buffer);

  size_t dwi_axis = 3;
  while (dwi_buffer.dim (dwi_axis) < 2) ++dwi_axis;
  INFO ("assuming DW images are stored along axis " + str (dwi_axis));

  Math::Matrix<cost_value_type> bmatrix;
  DWI::grad2bmatrix (bmatrix, grad);

  Math::Matrix<cost_value_type> binv (bmatrix.columns(), bmatrix.rows());
  Math::pinv (binv, bmatrix);

  int method = 1;
  Options opt = get_options ("method");
  if (opt.size()) method = opt[0][0];

  opt = get_options ("regularisation");
  cost_value_type regularisation = 5000.0;
  if (opt.size()) regularisation = opt[0][0];

  opt = get_options ("mask");
  Ptr<MaskBufferType> mask_buffer;
  Ptr<MaskBufferType::voxel_type> mask_vox;
  if (opt.size()){
    mask_buffer = new MaskBufferType (opt[0][0]);
    Image::check_dimensions (*mask_buffer, dwi_buffer, 0, 3);
    mask_vox = new MaskBufferType::voxel_type (*mask_buffer);
  }


  Image::Header dt_header (dwi_buffer);
  dt_header.set_ndim (4);
  dt_header.dim (3) = 6;
  dt_header.datatype() = DataType::Float32;
  dt_header.DW_scheme() = grad;

  OutputBufferType dt_buffer (argument[1], dt_header);

  InputBufferType::voxel_type dwi_vox (dwi_buffer);
  OutputBufferType::voxel_type dt_vox (dt_buffer);

  Image::ThreadedLoop loop ("estimating tensor components...", dwi_vox, 1, 0, 3);
  Processor processor (dwi_vox, dt_vox, mask_vox, bmatrix, binv, method, regularisation, loop.inner_axes()[0], dwi_axis);

  loop.run_outer (processor);
}

void verify_matrix (Math::Matrix<float>& in, const node_t num_nodes)
{
  if (in.rows() != in.columns())
    throw Exception ("Connectome matrix is not square (" + str(in.rows()) + " x " + str(in.columns()) + ")");
  if (in.rows() != num_nodes)
    throw Exception ("Connectome matrix contains " + str(in.rows()) + " nodes; expected " + str(num_nodes));

  for (node_t row = 0; row != num_nodes; ++row) {
    for (node_t column = row+1; column != num_nodes; ++column) {

      const float lower_value = in (column, row);
      const float upper_value = in (row, column);

      if (upper_value && lower_value && (upper_value != lower_value))
        throw Exception ("Connectome matrix is not symmetrical");

      if (!upper_value && lower_value)
        in (row, column) = lower_value;

      in (column, row) = 0.0f;

  } }
}

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 ()
{
    Image::Buffer<value_type> SH_data (argument[0]);
    Math::SH::check (SH_data);

    Options opt = get_options ("mask");

    std::unique_ptr<Image::Buffer<bool> > mask_data;
    if (opt.size())
        mask_data.reset (new Image::Buffer<bool> (opt[0][0]));

    opt = get_options ("seeds");
    Math::Matrix<value_type> dirs;
    if (opt.size())
        dirs.load (opt[0][0]);
    else {
        dirs.allocate (60,2);
        dirs = Math::Matrix<value_type> (default_directions, 60, 2);
    }
    if (dirs.columns() != 2)
        throw Exception ("expecting 2 columns for search directions matrix");

    opt = get_options ("num");
    int npeaks = opt.size() ? opt[0][0] : 3;

    opt = get_options ("direction");
    std::vector<Direction> true_peaks;
    for (size_t n = 0; n < opt.size(); ++n) {
        Direction p (Math::pi*to<float> (opt[n][0]) /180.0, Math::pi*float (opt[n][1]) /180.0);
        true_peaks.push_back (p);
    }
    if (true_peaks.size())
        npeaks = true_peaks.size();

    opt = get_options ("threshold");
    value_type threshold = -INFINITY;
    if (opt.size())
        threshold = opt[0][0];

    Image::Header header (SH_data);
    header.datatype() = DataType::Float32;

    opt = get_options ("peaks");
    std::unique_ptr<Image::Buffer<value_type> > ipeaks_data;
    if (opt.size()) {
        if (true_peaks.size())
            throw Exception ("you can't specify both a peaks file and orientations to be estimated at the same time");
        if (opt.size())
            ipeaks_data.reset (new Image::Buffer<value_type> (opt[0][0]));

        Image::check_dimensions (header, *ipeaks_data, 0, 3);
        npeaks = ipeaks_data->dim (3) / 3;
    }
    header.dim(3) = 3 * npeaks;
    Image::Buffer<value_type> peaks_data (argument[1], header);

    DataLoader loader (SH_data, mask_data.get());
    Processor processor (peaks_data, dirs, Math::SH::LforN (SH_data.dim (3)),
                         npeaks, true_peaks, threshold, ipeaks_data.get());

    Thread::run_queue (loader, Item(), Thread::multi (processor));
}

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();
    }

  }


}