C++ Vector::size方法代码示例

void run ()
{
  Image::Buffer<value_type> dwi_data (argument[0]);

  if (dwi_data.ndim() != 4)
    throw Exception ("dwi image should contain 4 dimensions");

  Math::Matrix<value_type> grad = DWI::get_valid_DW_scheme<value_type> (dwi_data);
  DWI::Shells shells (grad);
  // Keep the b=0 shell (may be used for normalisation), but force single non-zero shell
  shells.select_shells (true, true);

  Math::Matrix<value_type> DW_dirs = DWI::gen_direction_matrix (grad, shells.largest().get_volumes());

  Options opt = get_options ("lmax");
  lmax = opt.size() ? opt[0][0] : Math::SH::LforN (shells.largest().count());
  INFO ("calculating even spherical harmonic components up to order " + str (lmax));

  Math::Matrix<value_type> HR_dirs;
  Math::Matrix<value_type> HR_SHT;
  opt = get_options ("normalise");
  if (opt.size()) {
    normalise = true;
    opt = get_options ("directions");
    if (opt.size())
      HR_dirs.load (opt[0][0]);
    else
      DWI::Directions::electrostatic_repulsion_300 (HR_dirs);
    Math::SH::init_transform (HR_SHT, HR_dirs, lmax);
  }

  // set Lmax
  int i;
  for (i = 0; Math::SH::NforL(i) < shells.largest().count(); i += 2);
  i -= 2;
  if (lmax > i) {
    WARN ("not enough data for SH order " + str(lmax) + ", falling back to " + str(i));
    lmax = i;
  }
  INFO("setting maximum even spherical harmonic order to " + str(lmax));

  // Setup response function
  int num_RH = (lmax + 2)/2;
  Math::Vector<value_type> sigs(num_RH);
  std::vector<value_type> AL (lmax+1);
  Math::Legendre::Plm_sph<value_type>(&AL[0], lmax, 0, 0);
  for (int l = 0; l <= lmax; l += 2) sigs[l/2] = AL[l];
  Math::Vector<value_type> response(num_RH);
  Math::SH::SH2RH(response, sigs);

  opt = get_options ("filter");
  Math::Vector<value_type> filter;
  if (opt.size()) {
    filter.load (opt[0][0]);
    if (filter.size() <= response.size())
      throw Exception ("not enough filter coefficients supplied for lmax" + str(lmax));
    for (int i = 0; i <= lmax/2; i++) response[i] *= filter[i];
    INFO ("using initial filter coefficients: " + str (filter));
  }

  Math::SH::Transform<value_type> FRT_SHT(DW_dirs, lmax);
  FRT_SHT.set_filter(response);

  Image::Header qbi_header (dwi_data);
  qbi_header.dim(3) = Math::SH::NforL (lmax);
  qbi_header.datatype() = DataType::Float32;
  Image::Stride::set (qbi_header, Image::Stride::contiguous_along_axis (3, dwi_data));
  Image::Buffer<value_type> qbi_data (argument[1], qbi_header);


  opt = get_options ("mask");
  if (opt.size()) {
    Image::Buffer<bool> mask_data (opt[0][0]);
    Image::ThreadedLoop ("estimating dODFs using Q-ball imaging...", dwi_data, 0, 3)
      .run (DWI2QBI (FRT_SHT.mat_A2SH(), HR_SHT, shells), mask_data.voxel(), dwi_data.voxel(), qbi_data.voxel());
  }
  else 
    Image::ThreadedLoop ("estimating dODFs using Q-ball imaging...", dwi_data, 0, 3)
      .run (DWI2QBI (FRT_SHT.mat_A2SH(), HR_SHT, shells), dwi_data.voxel(), qbi_data.voxel());
}

void run () 
{

  Image::Header H (argument[0]);

  Image::Info info (H);
  info.set_ndim (3);
  Image::BufferScratch<bool> mask (info);
  auto v_mask = mask.voxel();

  std::string mask_path;
  Options opt = get_options ("mask");
  if (opt.size()) {
    mask_path = std::string(opt[0][0]);
    Image::Buffer<bool> in (mask_path);
    if (!Image::dimensions_match (H, in, 0, 3))
      throw Exception ("Input mask image does not match DWI");
    if (!(in.ndim() == 3 || (in.ndim() == 4 && in.dim(3) == 1)))
      throw Exception ("Input mask image must be a 3D image");
    auto v_in = in.voxel();
    Image::copy (v_in, v_mask, 0, 3);
  } else {
    for (auto l = Image::LoopInOrder (v_mask) (v_mask); l; ++l) 
      v_mask.value() = true;
  }

  DWI::CSDeconv<float>::Shared shared (H);

  const size_t lmax = DWI::lmax_for_directions (shared.DW_dirs);
  if (lmax < 4)
    throw Exception ("Cannot run dwi2response with lmax less than 4");
  shared.lmax = lmax;

  Image::BufferPreload<float> dwi (H, Image::Stride::contiguous_along_axis (3));
  DWI::Directions::Set directions (1281);

  Math::Vector<float> response (lmax/2+1);
  response.zero();

  {
    // Initialise response function
    // Use lmax = 2, get the DWI intensity mean and standard deviation within the mask and
    //   use these as the first two coefficients
    auto v_dwi = dwi.voxel();
    double sum = 0.0, sq_sum = 0.0;
    size_t count = 0;
    Image::LoopInOrder loop (dwi, "initialising response function... ", 0, 3);
    for (auto l = loop (v_dwi, v_mask); l; ++l) {
      if (v_mask.value()) {
        for (size_t volume_index = 0; volume_index != shared.dwis.size(); ++volume_index) {
          v_dwi[3] = shared.dwis[volume_index];
          const float value = v_dwi.value();
          sum += value;
          sq_sum += Math::pow2 (value);
          ++count;
        }
      }
    }
    response[0] = sum / double (count);
    response[1] = - 0.5 * std::sqrt ((sq_sum / double(count)) - Math::pow2 (response[0]));
    // Account for scaling in SH basis
    response *= std::sqrt (4.0 * Math::pi);
  }
  INFO ("Initial response function is [" + str(response, 2) + "]");

  // Algorithm termination options
  opt = get_options ("max_iters");
  const size_t max_iters = opt.size() ? int(opt[0][0]) : DWI2RESPONSE_DEFAULT_MAX_ITERS;
  opt = get_options ("max_change");
  const float max_change = 0.01 * (opt.size() ? float(opt[0][0]) : DWI2RESPONSE_DEFAULT_MAX_CHANGE);

  // Should all voxels (potentially within a user-specified mask) be tested at every iteration?
  opt = get_options ("test_all");
  const bool reset_mask = opt.size();

  // Single-fibre voxel selection options
  opt = get_options ("volume_ratio");
  const float volume_ratio = opt.size() ? float(opt[0][0]) : DWI2RESPONSE_DEFAULT_VOLUME_RATIO;
  opt = get_options ("dispersion_multiplier");
  const float dispersion_multiplier = opt.size() ? float(opt[0][0]) : DWI2RESPONSE_DEFAULT_DISPERSION_MULTIPLIER;
  opt = get_options ("integral_multiplier");
  const float integral_multiplier = opt.size() ? float(opt[0][0]) : DWI2RESPONSE_DEFAULT_INTEGRAL_STDEV_MULTIPLIER;

  SFThresholds thresholds (volume_ratio); // Only threshold the lobe volume ratio for now; other two are not yet used

  size_t total_iter = 0;
  bool first_pass = true;
  size_t prev_sf_count = 0;
  {
    bool iterate = true;
    size_t iter = 0;
    ProgressBar progress ("optimising response function... ");
    do {

      ++iter;

      {
        MR::LogLevelLatch latch (0);
        shared.set_response (response);
        shared.init();
//.........这里部分代码省略.........