C++ Header::DW_scheme方法代码示例

void run () {
  Point<> p[3];
  for (size_t arg = 0; arg < 3; ++arg) {
    std::vector<float> V = argument[arg]; 
    if (V.size() != 3) 
      throw Exception ("coordinates must contain 3 elements");
    p[arg][0] = V[0];
    p[arg][1] = V[1];
    p[arg][2] = V[2];
  }

  Image::Header header (argument[3]);
  float vox = Math::pow (header.vox(0)*header.vox(1)*header.vox(2), 1.0f/3.0f);
  Options opt = get_options ("vox");
  if (opt.size()) 
    vox = opt[0][0];


  Point<> d1 = p[1] - p[0];
  Point<> d2 = p[2] - p[0];
  d1.normalise();
  d2.normalise();

  d2 = d2 - d1 * d1.dot(d2);
  d2.normalise();

  Point<> d3 = d1.cross(d2);


  float min1 = std::numeric_limits<float>::infinity();
  float min2 = std::numeric_limits<float>::infinity();
  float max1 = -std::numeric_limits<float>::infinity();
  float max2 = -std::numeric_limits<float>::infinity();

  update_bounds (min1, min2, max1, max2, header, Point<int>(0,0,0), 0, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(0,0,0), 1, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(0,0,0), 2, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(1,0,0), 1, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(1,0,0), 2, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(0,1,0), 0, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(0,1,0), 2, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(0,0,1), 0, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(0,0,1), 1, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(1,1,0), 2, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(1,0,1), 1, p[0], d1, d2);
  update_bounds (min1, min2, max1, max2, header, Point<int>(0,1,1), 0, p[0], d1, d2);


  Point<> translation (p[0] + min1*d1 + min2*d2);

  header.transform()(0,0) = d1[0];
  header.transform()(1,0) = d1[1];
  header.transform()(2,0) = d1[2];

  header.transform()(0,1) = d2[0];
  header.transform()(1,1) = d2[1];
  header.transform()(2,1) = d2[2];

  header.transform()(0,2) = d3[0];
  header.transform()(1,2) = d3[1];
  header.transform()(2,2) = d3[2];

  header.transform()(0,3) = translation[0];
  header.transform()(1,3) = translation[1];
  header.transform()(2,3) = translation[2];

  header.set_ndim (3);
  header.dim (0) = (max1-min1) / vox;
  header.dim (1) = (max2-min2) / vox;
  header.dim (2) = 1;

  header.vox(0) = header.vox(1) = header.vox(2) = vox;

  header.datatype() = DataType::Bit;
  header.DW_scheme().clear();

  Image::Buffer<bool> roi_buffer (argument[4], header);
  Image::Buffer<bool>::voxel_type roi (roi_buffer);
  Image::LoopInOrder loop (roi);
  for (loop.start (roi); loop.ok(); loop.next (roi))
    roi.value() = true;

}

void run ()
{
  if (get_options ("norealign").size())
    Image::Header::do_not_realign_transform = true;

  const bool format     = get_options("format")        .size();
  const bool ndim       = get_options("ndim")          .size();
  const bool dimensions = get_options("dimensions")    .size();
  const bool vox        = get_options("vox")           .size();
  const bool dt_long    = get_options("datatype_long") .size();
  const bool dt_short   = get_options("datatype_short").size();
  const bool stride     = get_options("stride")        .size();
  const bool offset     = get_options("offset")        .size();
  const bool multiplier = get_options("multiplier")    .size();
  const bool comments   = get_options("comments")      .size();
  const bool properties = get_options("properties")    .size();
  const bool transform  = get_options("transform")     .size();
  const bool dwgrad     = get_options("dwgrad")        .size();

  Options opt = get_options ("export_grad_mrtrix");
  const std::string dw_out_mrtrix = opt.size() ? opt[0][0] : std::string();
  opt = get_options ("export_grad_fsl");
  std::string dw_out_fsl_bvecs, dw_out_fsl_bvals;
  if (opt.size()) {
    dw_out_fsl_bvecs = str(opt[0][0]);
    dw_out_fsl_bvals = str(opt[0][1]);
  }

  if ((dw_out_mrtrix.size() || dw_out_fsl_bvecs.size()) && (argument.size() > 1))
    throw Exception ("Can only export gradient table information to file if a single input image is provided");

  const bool print_full_header = (!(format || ndim || dimensions || vox || dt_long || dt_short || stride
                                  || offset || multiplier || comments || properties || transform || dwgrad)
                                  && dw_out_mrtrix.empty() && dw_out_fsl_bvecs.empty());


  for (size_t i = 0; i < argument.size(); ++i) {
    Image::Header header (argument[i]);

    if (format)     std::cout << header.format() << "\n";
    if (ndim)       std::cout << header.ndim() << "\n";
    if (dimensions) print_dimensions (header);
    if (vox)        print_vox (header);
    if (dt_long)    std::cout << (header.datatype().description() ? header.datatype().description() : "invalid") << "\n";
    if (dt_short)   std::cout << (header.datatype().specifier() ? header.datatype().specifier() : "invalid") << "\n";
    if (stride)     print_strides (header);
    if (offset)     std::cout << header.intensity_offset() << "\n";
    if (multiplier) std::cout << header.intensity_scale() << "\n";
    if (comments)   print_comments (header);
    if (properties) print_properties (header);
    if (transform)  std::cout << header.transform();
    if (dwgrad)     std::cout << header.DW_scheme();

    if (dw_out_mrtrix.size()) {
      if (!header.DW_scheme().is_set())
        throw Exception ("no gradient information found within image \"" + header.name() + "\"");
      header.DW_scheme().save (dw_out_mrtrix);
    }

    if (dw_out_fsl_bvecs.size()) {
      if (!header.DW_scheme().is_set())
        throw Exception ("no gradient information found within image \"" + header.name() + "\"");
      DWI::save_bvecs_bvals (header, dw_out_fsl_bvecs, dw_out_fsl_bvals);
    }

    if (print_full_header)
      std::cout << header.description();
  }

}

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

  }


}